7 views (last 30 days)

Show older comments

Hi everyone,

I have a set of experimental data of velocity fluctuation. The velocity plot involves two frequencies at 0.289 Hz and 0.578 Hz (See Figure 1 below). I would like to seprate the velocity into two waves at each frequency (0.289 Hz and 0.578 Hz) and to calculate their root mean square (rms) values. How can I do that? I tried with bandpass filter but the results show uniform distribution of the filtered velocity. Then, the rms can be easily calculated. Is this approach incorect? (See Figure 2 below)

Thanks a lot for your help.

Mathieu NOE
on 19 May 2021

hello again

I expanded a bit your code to try other bandpass filtering methods and also compare with my little fft code.

the fft peaks give a slightly lower amplitude compared to time signal rms computation (after bandpass filtering) , which is normal as it does not take the area within the same frequency range as the bandpass filter. Here there could be some further refinement if we wanted to proove Parseval's equation;

BTW your code and mine for bandpass filtering gives quite similar outputs. I overlaid the two outputs.

So nothing special to declare , the code is working fine IMHO

clearvars

clc

Fs = 20;

time_step = 1/Fs; %(s)

F = readtable('sl-5.csv', 'HeaderLines',1);

T = F(:,1);

L = table2array(T);

t = L.*time_step + time_step;

A=F{:,2};

A(isnan(A)) =0;

%%%%%%%% first frequency %%%%%%%%%%%%%

%% method 1

a = 0.1; % relative bandwith BP filter

low_1 = 0.289*(1-a/2);

high_1 = 0.289*(1+a/2);

[vel_filtered_1] = bandpass(A,[low_1 high_1],Fs);

%% method 2

a = 0.2; % relative bandwith BP filter

low_1 = 0.289*(1-a/2);

high_1 = 0.289*(1+a/2);

N = 4;

[numd,dend] = butter(N,[low_1 high_1]*2/Fs);

vel_filtered_1b = filtfilt(numd,dend,A);

% result_1 = [t vel_filtered_1];

%%%%%%%% second frequency %%%%%%%%%%%%%

%% method 1

a = 0.1; % relative bandwith BP filter

low_2 = 0.578*(1-a/2);

high_2 = 0.578*(1+a/2);

[vel_filtered_2] = bandpass(A,[low_2 high_2],Fs);

% result_2 = [t vel_filtered_2];

%% method 2

a = 0.2; % relative bandwith BP filter

low_2 = 0.578*(1-a/2);

high_2 = 0.578*(1+a/2);

N = 4;

[numd,dend] = butter(N,[low_2 high_2]*2/Fs);

vel_filtered_2b = filtfilt(numd,dend,A);

figure(1)

subplot(211)

plot(t, vel_filtered_1,t, vel_filtered_1b)

h= gca;

set(h,'YMinorTick','off','XMinorTick','on','LineWidth',2,'Box','on','FontSize',40);

ylabel('v (m/s)','FontSize',10);

xlabel('Time (s)','FontSize',10);

subplot(212)

plot(t, vel_filtered_2,t, vel_filtered_2b)

h= gca;

set(h,'YMinorTick','off','XMinorTick','on','LineWidth',2,'Box','on','FontSize',40);

ylabel('v (m/s)','FontSize',10);

xlabel('Time (s)','FontSize',10);

% generate rms values

vel_filtered_1_rms = sqrt(mean(vel_filtered_1.^2));

vel_filtered_2_rms = sqrt(mean(vel_filtered_2.^2));

disp(['Velocity at first frequency (0.289 Hz) : ' num2str(vel_filtered_1_rms) ' m/s rms']);

disp(['Velocity at second frequency (0.578 Hz) : ' num2str(vel_filtered_2_rms) ' m/s rms']);

% comparison with averaged windowed fft computation :

nfft = round(length(A)/8);

Overlap = 0.75;

[freq_vector,fft_spectrum1] = myfft_peak(A, Fs, nfft, Overlap);

% convert peak to rms values (divide by sqrt(2))

fft_spectrum1 = fft_spectrum1/1.414;

% zoom in 0 - 1 Hz frequency range

ind = find(freq_vector>0.1 & freq_vector<0.8);

freq_vector = freq_vector(ind);

fft_spectrum1 = fft_spectrum1(ind);

[val,locs] = findpeaks(fft_spectrum1,'NPeaks',2,'SortStr','descend');

freqpeaks = freq_vector(locs)

val

figure(2),plot(freq_vector,fft_spectrum1,'b',freqpeaks,val,'rd');

ylabel('v (m/s)','FontSize',10);

xlabel('Frequency (Hz)','FontSize',10);

% xlim([0.1 0.8]);

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

function [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)

% FFT peak spectrum of signal (example sinus amplitude 1 = 0 dB after fft).

% Linear averaging

% signal - input signal,

% Fs - Sampling frequency (Hz).

% nfft - FFT window size

% Overlap - buffer percentage of overlap % (between 0 and 0.95)

[samples,channels] = size(signal);

% fill signal with zeros if its length is lower than nfft

if samples<nfft

s_tmp = zeros(nfft,channels);

s_tmp((1:samples),:) = signal;

signal = s_tmp;

samples = nfft;

end

% window : hanning

window = hanning(nfft);

window = window(:);

% compute fft with overlap

offset = fix((1-Overlap)*nfft);

spectnum = 1+ fix((samples-nfft)/offset); % Number of windows

% % for info is equivalent to :

% noverlap = Overlap*nfft;

% spectnum = fix((samples-noverlap)/(nfft-noverlap)); % Number of windows

% main loop

fft_spectrum = 0;

for i=1:spectnum

start = (i-1)*offset;

sw = signal((1+start):(start+nfft),:).*(window*ones(1,channels));

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 % 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;

end

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

Start Hunting!