fft finding oscillation frequency and strouhal number for lift
26 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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
0 commentaires
Réponse acceptée
Mathieu NOE
le 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 commentaires
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Spectral Measurements dans Help Center et File Exchange
Produits
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!