28 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

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Mathieu NOE
on 29 Nov 2020

hello Megi

below I have a bit modified your code to improve the fft computation

and also compute the frequency of the first peak in the fft (plotted as a red star)

the frequency of the peak is around 7 Hz, not 1.9531 Hz

the value you mentionned (1.9531 Hz) is my fft frequency resolution, which is something different;

the fft spectrum is computed on a frequency vector and its resolution is related to fs and nfft by df = fs / nfft

that means the frequency increment is df = fs / nfft

in your case the nfft is equal to the length of the record, so we cannot reduce (refine) the frequency resolution;

you should (if needed) increase the time duration of your record

you could also reduce the sampling frequency because 2000 Hz is very fast for analysing signals below 10 Hz

you could easily reduce fs to 100 Hz, that would reduce the record file size and post processing time;

a=xlsread('lift.xlsx');

x=a(:,1);

y=a(:,2);

N=length(a); % sample number

fs=2000; % [Hz] sampling frequency

nfft = N; % no averaging : the fft length = the samples qty

% X=2*fft(y)/N;

window = hanning(nfft);

% fft_spectrum=abs(fft(y.*window))*4/nfft; % the fft amplitude is correct in this fomula if you stick to the hanning window

% one sidded fft spectrum % Select first half

if rem(nfft,2) % nfft odd

select = (1:(nfft+1)/2)';

else

select = (1:nfft/2+1)';

end

fft_spectrum = fft_spectrum(select);

freq_vector = (select - 1)*fs/nfft;

% let's find frequency of first peak

[peaks,loc] = findpeaks(fft_spectrum);

freq_of_first_peak = freq_vector(loc(1));

semilogx(freq_vector,fft_spectrum,'b',freq_of_first_peak,peaks(1),'*r');grid

xlabel('frequency');

ylabel('fft magnitude');

title('lift');

% strouhal number

St = freq_of_first_peak*D/U;

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

Start Hunting!
## 0 Comments

Sign in to comment.