b0 + b1*z^-1 + b2*z^-2 H(z) = _________________________ 1 + a1*z^-1 + a2*z^-2
Shock Response Spectrum (SRS) analysis was developed as a standard data processing method in the early 1960’s. Firstly SRS was used by U.S. Department Of Defense. Now this signal processing method is standardized by ISO 18431-4 Mechanical vibration and shock — Signal processing — Part 4: Shock response spectrum analysis
Example: system should take sweep sine, and simulate forced response
#ifndef SDOF_DIGITAL_H #define SDOF_DIGITAL_H #include &l tiostream > #include &l cmath > #include &l string > #include &l sstream > //using namespace std; /* DIGITAL H(w) = FRF(f) ===== > to H(z) _____ _____ ____ ______ * . Q = 1/2D / ____| __ \ / __ \| ____| | / \ | (___ | | | | | | | |__ |..../ | \ \___ \| | | | | | | __| | fr \............... [N] ____) | |__| | |__| | | |-------------------------> |_____/|_____/ \____/|_| || fmin fmax DIGITAL SYSTEM ::: sdof_digital.h // Difference equation for computation of extreme response spectrum // Explanation *ISO 18431-4 Standard Response of a linear Single Degree of Freedom System (SDOF system) to vibration, according to its natural frequency, for a given damping ratio. The response is described here by the relative movement of the mass of this system in relation to its support. //x2 = β0 x1 k + β1 x1(k −1)+ β2 x1(k − 2) − α1x2(k −1) − α2 x2(k − 2) // RECALCULATING SDOF filter coefficients /* Shock Response Spectrum is defined as the response to a given acceleration acting at a set of mass-damper-spring oscillators wnS / Q + wn^2 A2/A1 = ______________________ --> transformed to Z domain .... s^2 + (wn/Q)S +1 b0 + b1*z^-1 + b2*z^-2 H(z) = _________________________ 1 + a1*z^-1 + a2*z^-2 Jan Gusic, 7.5, 2021 */
/* DIGITAL SYSTEM FORM *************************************************************** Sdof_digital is a calss for calculating coefficients, which are in turn returned to the member structure C ( coefficients) ... this structure holds the difference equation coeffs. C --- | a1 | a2 | b0 | b1 | b2 Now returning the structure to the public structure of 'deq' ( diff. equations interface) coefficients are taken by the internal function deq MyEquation; Sdof_digital SDOF; which returns structure C after calling returnCoefficients() ( Ex. for first coefficients can be called SDOF.returnCoefficients().a1, .. ) Now, put the SDOF to Diff equation for single value as... MyEquation.process(Sdof_digital::C& Coeff, 0 ) And check the coefficients passed by reference from Sdof_digital instance as MyEquation.type() cout << SDOF.type(); //they should be the same.. FOR EACH signal value in signal_in MyEquation.process(SDOF.returnCoefficients(), signal_in(i)) END */
class Sdof_digital { private: double PI = 3.14159265358979323846; //| pi double fs; //| sampling f double fr; //| resonance double Q; //| Q double Ts; //| samp. time double wn; //| nat, freq. // BiQuad filter coefficients double A, B; double b0,b1,b2; double a1, a2; public: struct C{ double b0,b1,b2, a1, a2; }; Sdof_digital(/* default == Bypass */); //| Bypass sinal Sdof_digital(double, double, double); //| Digital system initialization ~Sdof_digital(); std::string type(); //| write coefficients to the screen C returnCoefficients(); //| function returns C structure after reading them .. };
Sdof_digital::Sdof_digital() { //bypass coefficients for xin = xout .. no processing // with this fitler * off b0 = 1; b1 = b2 = 0; a1 = a2 = 0; } std::string Sdof_digital::type() { std::stringstream sout; sout << "| fs:" << this->fs << "| fr:" << fr << "| Q: " << Q << "\n"; sout << "| a0:" << this->a1 << " | a1:" << this->a2 << "\n"; sout << "| b0:" << this->b0 << "| b1:" << this->b1 << "| b2:"<< this->b2 <<"\n"; return sout.str(); } Sdof_digital::~Sdof_digital() { } Sdof_digital::Sdof_digital(double Fs, double fr, double Q) { //| sampling //| | frequency //| | Q this->fs = Fs; this->Ts = 1/fs; this->wn = 2*PI*fr; this->fr = fr; this->Q = Q; this->A = wn*Ts/2/Q; this->B = wn*Ts*sqrt(1-1/4/Q/Q); this->b0 = 1-exp(-A)*sin(B)/B; this->b1 = 2*exp(-A)*(sin(B)/B-cos(B)); this->b2 = exp((-2)*A)-exp(-A)*sin(B)/B; this->a1 = (-2)*exp((-1)*A)*cos(B); this->a2 = exp((-2)*A); } Sdof_digital::C Sdof_digital::returnCoefficients() { C coeff; coeff.b0 = this->b0; coeff.b1 = this->b1; coeff.b2 = this->b2; coeff.a1 = this->a1; coeff.a2 = this->a2; return coeff; } // Difference equation struct deq { // coefficients double b0, b1, b2; double a1,a2; //input history double x1,x2; //output history double y1, y2; // procesurianje ::: Biquad filter normalized double process(const Sdof_digital::C& COEFF,double x0){ b0 = COEFF.b0; b1 = COEFF.b1; b2 = COEFF.b2; a1 = COEFF.a1; a2 = COEFF.a2; double y0 = b0*x0 + b1*x1 + b2*x2 - a1*y1 -a2*y2; x2 = x1; x1 = x0; y2 = y1; y1 = y0; return y0; } // type the filtering coefficients to the console .. void type() { std::cout << "\n\n Difference equation setup:::"; std::cout <<"| bo:" << b0 << "| b1: " << b1 << "|b2: " << b2 << std::endl; std::cout <<"| a1:" << a1 << "| a2:" << a2 << "... " << std::endl; } }; #endif
#include "sdof_digital.h" #include "signal_gen.h" using namespace std; // Functionality of sdof_digital.h // herewith confrimed .. // // int main() { //(1) // Define Digital system :: Sdof_digital MySdof(fs,fr,Q) Sdof_digital sys1; Sdof_digital sys2(21000, 400, 40); cout << sys2.type(); //(2) // Define Digital process :: deq CH_X; i.e. Turn the system to H(z) , get coefficients // then CH_X(const Sdof_digital::C&, signal_input ) deq CH1; CH1.process(sys2.returnCoefficients(),1); CH1.type(); //(3) // Define input signal from signal_gen.h signal_gen sineWave; sineWave.sweep_linear(21000,1, 1e3, 1, 0,1); vector <double> x = sineWave.getSignal(); //(4) // define the output vector to receive processed series ... vector < double > y; double yy; //(5) // process the series // x2 = β0 x1 k + β1 x1(k −1)+ β2 x1(k − 2) − α1x2(k −1) − α2 x2(k − 2) for (size_t i = 0; i < x.size(); i++) { yy = CH1.process(sys2.returnCoefficients(), x[i]); y.push_back(yy); } //(6) write to file for Gnuplot.. ofstream results("mySys_test.dat"); for(size_t i = 0; i < y.size(); i++){ results << i << "\t" << x[i] << "\t" << y[i] <<endl; } results.close(); return 0; }
(1) Input sinewave (2) input and forced response (3) Zooom into max. forced response region
Forced response and input sweep sinewave