fft finding oscillation frequency and strouhal number for lift

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

Answer 1

[Mathieu NOE]

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

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

[Megi Clen]

Hello Mathieu,

Thanks for the answer. Your script seems very complicated to me as a beginner. I tried to copy and run it, but there is an error coming from the line below. Sorry but what could be the reason? when using your script, do I need to change anything particularly or I can use all of it?

Note: Also, sorry for the script picture, it is my first time to ask question here. Next time, I will be careful.

Regards,

Megi

A=readmatrix('lift.xlsx');

[Mathieu NOE]

hello Megi

there should be no modification needed to run the code

sure , it takes time to master matlab (even myself have still a lot to learn)... but don't be afraid of a complicated code, you'll find all needed explanations at the right moment.

so which release of matlab do you use ?

there is an older function to read xls file which is readxls

if you have a recent matlab, you should have also readtable

can you check please ?

[Mathieu NOE]

sorry

I meant xlsread and not readxls

so this code should also work

% A = readmatrix('lift.xlsx');
A = xlsread('lift.xlsx');

[Megi Clen]

Hi Mathieu,

Thanks. I am using MATLAB R2018a version. I tried the code again but still didn't go nowhere. This time, the programme is giving following error:

error using plot

vectors must be the same length.

error in grid (line 14)

plot(f, abs (z))

error in matlabInfo (line 53)

figure (1), plot (freq, sensor_spectrum_dB, 'b'); grid

Do you have any idea what is wrong about it?

Thanks,

Megi

[Mathieu NOE]

hello

sorry for that

I made a few modifications in the my_fft function

hope it works now

m file is attached

please let me know

[Megi Clen]

Thank you Mathieu but new script didn't work as well. maybe I should just leave it. thanks for your time and help.

[Mathieu NOE]

hum

that sounds bizarre to me

could you copy the error message and also what are the dimensions of the data that cause the trouble ? (type who in the command window)

[Megi Clen]

Ok, the error is below:

error using plot

vectors must be the same length.

error in grid (line 14)

plot (f, abs(z))

error in simple_demo_ok (line 72)

figure (1), semilogx (freq, sensor_spectrum_dB, 'b'); grid

[Mathieu NOE]

hmmm

I wonder if you have another function called grid which is not the native matlab function and interfere with what I wanted to do;

because grid does not have a line 14 which does plot (f, abs(z))

so try again with grid removed from the code

FYI , when I edit grid , this is what it shows :

function grid(varargin)
%GRID   Grid lines.
%   GRID ON adds major grid lines to the current axes.
%   GRID OFF removes major and minor grid lines from the current axes. 
%   GRID MINOR toggles the minor grid lines of the current axes.
%   GRID, by itself, toggles the major grid lines of the current axes.
%   GRID(AX,...) uses axes AX instead of the current axes.
%
%   GRID sets the XGrid, YGrid, and ZGrid properties of
%   the current axes. If the axes is a polar axes then GRID sets
%   the ThetaGrid and RGrid properties. If the axes is a geoaxes, then GRID
%   sets the Grid property.
%
%   AX.XMinorGrid = 'on' turns on the minor grid.
%
%   See also TITLE, XLABEL, YLABEL, ZLABEL, AXES, PLOT, BOX, POLARAXES.

%   Copyright 1984-2019 The MathWorks, Inc.

% To ensure the correct current handle is taken in all situations.

isaxarr=matlab.graphics.chart.internal.objArrayDispatch(@grid,varargin{:});
if isaxarr
    return
end

narginchk(0,2)

opt_grid = 0;
if nargin == 0
    ax = gca;
    
    % Chart subclass support
    if isa(ax,'matlab.graphics.chart.Chart')
        grid(ax);
        return
    end
else
    arg1=varargin{1};

    if matlab.graphics.internal.isCharOrString(arg1)
        % update string input to OnOffSwitchState if applicable
        if nargin == 2
            error(message('MATLAB:grid:FirstArgAxes'))
        end
        ax = gca;
        
        if strcmpi(arg1, 'minor')
            opt_grid = 'minor';
        else
            try
                opt_grid = matlab.lang.OnOffSwitchState(lower(arg1));
            catch
                error(message('MATLAB:grid:UnknownOption'));
            end
        end
        
        % Chart subclass support
        if isa(ax,'matlab.graphics.chart.Chart')
            grid(ax,opt_grid);
            return
        end
    else
        if ~any(isgraphics(arg1,'axes') ...
                | isgraphics(arg1,'polaraxes') ...
                | isgraphics(arg1,'geoaxes'),'all')
            error(message('MATLAB:grid:FirstArgAxes'));
        end
        ax = arg1;
        
        % check for string option
        if nargin == 2
            in2 = varargin{2};
            if isscalar(in2) && isnumeric(in2) && in2 == 0
               opt_grid = 0;
            elseif strcmpi(in2, 'minor')
                opt_grid = 'minor';
            else
                try
                    opt_grid = matlab.lang.OnOffSwitchState(in2);
                catch
                    error(message('MATLAB:grid:UnknownOption'));
                end
            end
        end
    end
end

if isempty(opt_grid)
    error(message('MATLAB:grid:UnknownOption'));
end

names = get(ax,'DimensionNames');
xgrid = [names{1} 'Grid'];
ygrid = [names{2} 'Grid'];
zgrid = [names{3} 'Grid'];
xminorgrid = [names{1} 'MinorGrid'];
yminorgrid = [names{2} 'MinorGrid'];
zminorgrid = [names{3} 'MinorGrid'];

matlab.graphics.internal.markFigure(ax);

%---Check for bypass option
if isappdata(ax,'MWBYPASS_grid')
   mwbypass(ax,'MWBYPASS_grid',opt_grid);
elseif isgraphics(ax,'geoaxes')
    % geoaxes only has 1 Grid property and does not have a minor grid
    if isnumeric(opt_grid) % opt_grid == 0
        set(ax,'Grid', ~get(ax,'Grid'));
    elseif (strcmp(opt_grid, 'minor'))
        error(message('MATLAB:Chart:UnsupportedArgument','grid','geoaxes'));
    else % opt_grid == OnOffSwitchState on/off
        set(ax,'Grid', opt_grid)
    end
elseif strcmp(opt_grid, 'minor')
    set(ax,xminorgrid, ~get(ax,xminorgrid));
    set(ax,yminorgrid, ~get(ax,yminorgrid));

    if hasZProperties(handle(ax))
        set(ax,zminorgrid, ~get(ax,zminorgrid));
    end
elseif isequal(opt_grid, 0)&& ~isa(opt_grid, 'matlab.lang.OnOffSwitchState')
    set(ax,xgrid,~get(ax,xgrid));
    set(ax,ygrid,~get(ax,ygrid));

    if hasZProperties(handle(ax))
        set(ax,zgrid,~get(ax,zgrid));
    end
else
    set(ax, xgrid, opt_grid, ygrid, opt_grid);
    if ~opt_grid
        set(ax, xminorgrid, opt_grid, yminorgrid, opt_grid);
    end
    if hasZProperties(handle(ax))
        set(ax, zgrid, opt_grid);
        if ~opt_grid
            set(ax,zminorgrid, opt_grid);
        end
    end
end

[Megi Clen]

Ok, thanks. I will try again and let you know.

[Megi Clen]

Hi Mathieu, I want to let you know that your scripts are working fine. Sorry, the problem I was having because of a file named "grid.m" in my folder. After deleting it, I got the plot. Thank you so much.

Also, regarding Strouhal number (St) graph, can I use the frequency obtained from the graph (f=1.9531 Hz) for the script below? i am asking because when I use the sampling frequency as 2000 Hz, I am taking very high values of St as in the picture. However, I need St values smaller than 1 for my case. Do you have any idea for this?

Note: St= f*D/U (f=oscillation frequency, D=cylinder diameter, U=flow velocity)

Regards, Megi.

a=xlsread('lift.xlsx');
x=a(:,1);
y=a(:,2);
N=length(a); % sample number
fs=2000; % [Hz] sample frequency
X=2*fft(y)/N;
fl=linspace(0,fs,N);
plot(fl,abs(X));
xlabel('strouhal number');
ylabel('magnitude');
title('lift');

[Mathieu NOE]

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;

[Megi Clen]

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.

[Mathieu NOE]

you're welcome

good luck for the future