fft finding oscillation frequency and strouhal number for lift

26 views (last 30 days)
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

Accepted Answer

Mathieu NOE
Mathieu NOE on 25 Nov 2020
hello
it's difficult to test your code if you share it only as a picture
it looks ok but if you want to double check this is a code I regularly share for fft spectral analysis
looks like the oscillation frequency is around 8 Hz
nb : your time vector is not equally spaced in your file , maybe due to rounding somewhere
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 1024; %
OVERLAP = 0.85;
% spectrogram dB scale
spectrogram_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 data
A = readmatrix('lift.xlsx');
time = A(:,1); % NB : time stamps are not equally spaced
signal = 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 needed
if 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))');
end
figure(1),plot(freq,sensor_spectrum_dB,'b');grid
title(['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 needed
if 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 spectrogram
figure(2);
imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
axis('xy');colorbar('vert');grid
title([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(:));
end
function [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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  14 Comments
Megi Clen
Megi Clen on 29 Nov 2020
Hi Mathieu. I got the new plot as shown in the picture. Also, I must tell you that the new script needs a little modification, when % symbol in the line 7 and 9 is removed, it is working fine.
Many thanks for your help and information.
Have a nice day, Megi.

Sign in to comment.

More Answers (0)

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!