# MATLAB: Fft finding oscillation frequency and strouhal number for lift

fftMATLABoscillation frequency and strouhal number

Hello all,
I have a data which comes from an oscillating plate in a laminar flow. (please see the attachment). my sampling frequency is 2000 Hz. I am trying find out frequency magnitude and strouhal number graphs with fft. For the frequency plot, I would like to neglect all frequencies above the Nyquist frequency. I have looked at the examples on the page, followed them and created my scripts as below but at some points I became confused and not sure about my script is right or not. So, I am wondering does any of you help me to solve this problem ?
Thanks,
Megi

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FFT parameters%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%NFFT = 1024;    % OVERLAP = 0.85;% spectrogram dB scalespectrogram_dB_scale = 80;  % dB range scale (means , the lowest displayed level is XX dB below the max level)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% options %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if you are dealing with acoustics, you may wish to have A weighted% spectrums % option_w = 0 : linear spectrum (no weighting dB (L) )% option_w = 1 : A weighted spectrum (dB (A) )option_w = 0;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% load signal%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% test dataA = readmatrix('lift.xlsx');time = A(:,1); % NB : time stamps are not equally spacedsignal = A(:,2);dt = mean(diff(time));Fs = 1/dt;Fs = 2000; % had to force the value because time stamps are not equally spaced% % decimate% decim = 12;% signal = decimate(signal,decim);% Fs = Fs/decim;% %[data,Fs]=wavread('Approach_Gear_Drop_Aft Ctr.wav ');  %(older matlab)% % or % [data,Fs]=audioread('myWAVaudiofile.wav'); %(newer matlab)% channel = 1;% signal = data(:,channel);% samples = length(signal);dt = 1/Fs;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% display 1 : averaged FFT spectrum%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[freq, sensor_spectrum] = myfft_peak(signal,Fs,NFFT,OVERLAP);% convert to dB scale (ref = 1)sensor_spectrum_dB = 20*log10(sensor_spectrum);% apply A weigthing if neededif option_w == 1    pondA_dB = pondA_function(freq);    sensor_spectrum_dB = sensor_spectrum_dB+pondA_dB;    my_ylabel = ('Amplitude (dB (A))');else    my_ylabel = ('Amplitude (dB (L))');endfigure(1),plot(freq,sensor_spectrum_dB,'b');gridtitle(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);xlabel('Frequency (Hz)');ylabel(my_ylabel);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% display 2 : time / frequency analysis : spectrogram demo%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[sg,fsg,tsg] = specgram(signal,NFFT,Fs,hanning(NFFT),floor(NFFT*OVERLAP));  % FFT normalisation and conversion amplitude from linear to dB (peak)sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg));     % NB : X=fft(x.*hanning(N))*4/N; % hanning only% apply A weigthing if neededif option_w == 1    pondA_dB = pondA_function(fsg);    sg_dBpeak = sg_dBpeak+(pondA_dB*ones(1,size(sg_dBpeak,2)));    my_title = ('Spectrogram (dB (A))');else    my_title = ('Spectrogram (dB (L))');end% saturation of the dB range : min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;% plots spectrogramfigure(2);imagesc(tsg,fsg,sg_dBpeak);colormap('jet');axis('xy');colorbar('vert');gridtitle([my_title ' / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);xlabel('Time (s)');ylabel('Frequency (Hz)');%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%รนfunction pondA_dB = pondA_function(f)	% dB (A) weighting curve	n = ((12200^2*f.^4)./((f.^2+20.6^2).*(f.^2+12200^2).*sqrt(f.^2+107.7^2).*sqrt(f.^2+737.9^2)));	r = ((12200^2*1000.^4)./((1000.^2+20.6^2).*(1000.^2+12200^2).*sqrt(1000.^2+107.7^2).*sqrt(1000.^2+737.9^2))) * ones(size(f));	pondA = n./r;	pondA_dB = 20*log10(pondA(:));endfunction  [f,fft_spectrum] = myfft_peak(s, fs, nfft, Overlap)    %FFT peak spectrum of signal  (example sinus amplitude 1   = 0 dB after fft).    %   s - input signal,     %   fs - Sampling frequency (Hz).    %   spectsize - FFT window size    %   spectnum - number of windows to analyze    f = [0:fs/nfft:fs/2];    f = f(:);    dt = 1/fs;    samples = length(s);    % fill s with zeros if length is lower than nfft    if samples<nfft        s_tmp = zeros(nfft,1);        s_tmp((1:samples)) = s;        s = s_tmp;    end    % window : hanning    window = hanning(nfft);    window = window(:);    % compute fft with overlap     offset = floor((1-Overlap)*nfft);    spectnum = 1+ floor((length(s)-nfft)/offset); % nb of averages    fft_spectrum = 0;        % main loop        for i=1:spectnum            start = (i-1)*offset;            sw = s((1+start):(start+nfft)).*window;            fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft);     % X=fft(x.*hanning(N))*4/N; % hanning only         end        fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling    % one sidded fft spectrum (index= (1:L/2+1))    fft_spectrum = fft_spectrum(1:nfft/2+1);    end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%