​

   

    

     

      

     
     

      

       clear, close all
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      

       %% initialize parameters
      
     
    

     

      

     
     

      

       samplerate=
       500; % 
       in Hz
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      

       fstep=
       1;   % frequency 
       step 
       for wavelet
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      

       %% generate simulated signals 
       with 
       step changes 
       in frequency
      
     
    

     

      

     
     

      

       data = csvread(
       '3_1_link6_28_5_30min.csv'); % input the signal from the Excle
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      

       data = data
       '; % change the signal from column to row
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      

       N = length(data);   % calculate the length 
       of the data
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      

       taxis = [
       1:N]/samplerate;   % time axis 
       for whole data length
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      

       figure, 
      
     
    

     

      

     
     

      

       plot(taxis,data),xlim([taxis(
       1) taxis(
       end)])
      
     
    

     

      

     
     

      

       xlabel(
       'Time (s)')
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      

       %% Time-frequency analysis (CWT, morlet wavelet)
      
     
    

     

      

     
     

      

       spec = tfa_morlet(data, samplerate, 
       1, 
       250, fstep);
      
     
    

     

      

     
     

      

       faxis=[
       1:fstep:
       250];
      
     
    

     

      

     
     

      

       Mag=abs(spec);     % 
       get spectrum magnitude
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      

       im = figure(
       'color',[1 1 1]);
      
     
    

     

      

     
     

      

       imagesc(taxis,faxis,Mag)   % plot spectrogram 
       as an image
      
     
    

     

      

     
     

      

       colorbar
      
     
    

     

      

     
     

      

       axis([taxis(
       1) taxis(
       end) faxis(
       1) faxis(
       end)])
      
     
    

     

      

     
     

      

       xlabel(
       'Time (s)')
      
     
    

     

      

     
     

      

       ylabel(
       'Frequency (Hz)')
      
     
    

     

      

     
     

      

       title(
       'Time-frequency analysis (CWT)')
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      

       saveas(im,
       'CWT_1.bmp')
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      

       function TFmap = tfa_morlet(td, fs, fmin, fmax, fstep)
      
     
    

     

      

     
     

      

       TFmap = [];
      
     
    

     

      

     
     

      

       for fc=fmin:fstep:fmax
      
     
    

     

      

     
     

      

           MW = MorletWavelet(fc/fs);  % calculate the Morlet Wavelet 
       by giving the central freqency
      
     
    

     

      

     
     

      

           cr = conv(td, MW, 
       'same'); % convolution
      
     
    

     

      

     
     

      
    
      
     
    

     

      

     
     

      

           TFmap = [TFmap; abs(cr)];
      
     
    

     

      

     
     

      

       end
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      

       function MW = MorletWavelet(fc)
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      

       F_RATIO = 
       7;    % frequency ratio (number 
       of cycles): fc/sigma_f, should be greater than 
       5
      
     
    

     

      

     
     

      

       Zalpha2 = 
       3.3;  % value 
       of Z_alpha/
       2, 
       when alpha=
       0.001
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      

       sigma_f = fc/F_RATIO;
      
     
    

     

      

     
     

      

       sigma_t = 
       1/(
       2*pi*sigma_f);
      
     
    

     

      

     
     

      

       A = 
       1/sqrt(sigma_t*sqrt(pi));
      
     
    

     

      

     
     

      

       max_t = ceil(Zalpha2 * sigma_t);
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      

       t = -max_t:max_t;
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      

       %MW = A * exp((-t.^
       2)/(
       2*sigma_t^
       2)) .* exp(
       2i*pi*fc*t);
      
     
    

     

      

     
     

      

       v1 = 
       1/(-
       2*sigma_t^
       2);
      
     
    

     

      

     
     

      

       v2 = 
       2i*pi*fc;
      
     
    

     

      

     
     

      

       MW = A * exp(t.*(t.*v1+v2));