Search code examples
matlabsignal-processingfftodewavelet

FFT of time series generated from a nonlinear dynamical system


       close all;clc;
        global  a b c
         a=0.2;
         b=0.4;
         c=5.7;
     ts=0:.01:4000;
     z0=[1 0 1]; 
     opt=odeset('RelTol',10e-12);
     [t,z]= ode45('System', ts, z0,opt);
     Fs = 1000;   
     x=z(:,1);
     nfft = 2^nextpow2(length(x));
     Pxx = abs(fft(x,nfft)).^2/length(x)/Fs;
     Hpsd = dspdata.psd(Pxx(1:length(Pxx)/2),'Fs',Fs/10); 
     figure()
     plot(Hpsd)

    function zdot=System(t,z) 
      global a b c
   zdot=[ -(z(2)+z(3));
           z(1)+a*z(2);
           b+z(1)*z(3)-c*z(3)];

The code above is that of Rossler nonlinear (chaotic) dynamical system. The integration time step is 0.01 that is a parameter to ode45 solver. I am trying to plot the fft and the scalogram (wavelet transform). But I don't know how to set the the sampling and nyquist fequency for these types of systems. I just assumed sampling frequency to be Fs = 1000.....it could be 100 as well. This is what I am unsure of.

Can somebody please help in explaining

(1) what is the sampling frequency and nyquist frequency for these types of dynamical systems

(2) how to obtain the fft and the scalogram image from wavelet transform of z


Solution

  • Assuming that there's a continuous-time dynamical system, whose governing system of differential equations, represented in suitable system of units such as SI system, produces an output as its solution, say y(t) for simplicity, then you convert this continuous-time function y(t) into a discrete-time sequence y[n] by sampling it properly.

    The most typical sampling method is the uniform sampling: y[n] = y(tn) = y(n*Ts).

    where Ts is the sampling period in seconds. This means that the sample values y[n] of the discrete-time sequence is obtained by evaluating the continuous-time function y(t) at the time-stamps given by : tn = n*Ts .

    With this sampling method, the sampling frequency is Fs = 1/Ts in Hertz, and the associated Nyquist frequency Fn = Fs/2.

    In your code, the ODE45 implements a numerical procedure to solve the presented system of differential equations, evaluated at the time-stemps indicated by the vector: ts=0:.01:4000; from 0 to 4000 seconds (assuming second as the unit of time)

    This time-vector ts indicates the sampling of the continuous-time solution y(t), and the step-size of this vector is, what is defined above as, the sampling period Ts.

    Therefore, in your specific example the sampling period is 0.01 seconds and sampling frequency is 100 Hz and Nyquist frequency is 50 Hz.