# separate signals and calculate their rms at different frequencies

7 views (last 30 days)
Quang Duy Nguyen on 17 May 2021
Commented: Quang Duy Nguyen on 19 May 2021
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.
##### 2 CommentsShowHide 1 older comment
Quang Duy Nguyen on 18 May 2021
Hi Mathieu, thank you for your reply.
I used bandpass filters within +- 5% from the centre values (0.289 Hz and 0.578 Hz). I don't know why the filters return the results quite uniform and the time-average values of the filtered are zeros. Please find the code and the data below.
Thanks
clear all
clc
Fs = 20;
time_step = 1/Fs; %(s)
T = F(:,1);
L = table2array(T);
t = L.*time_step + time_step;
A=F{:,2};
A(isnan(A)) =0;
low_1 = 0.289*(1-0.05);
high_1 = 0.289*(1+0.05);
[vel_filtered_1] = bandpass(A,[low_1 high_1],Fs);
result_1 = [t vel_filtered_1];
low_2 = 0.578*(1-0.05);
high_2 = 0.578*(1+0.05);
[vel_filtered_2] = bandpass(A,[low_2 high_2],Fs);
result_2 = [t vel_filtered_2];
figure
subplot(211)
plot(t, vel_filtered_1)
h= gca;
set(h,'YMinorTick','off','XMinorTick','on','LineWidth',2,'Box','on','FontSize',40);
ylabel('v (m/s)','FontSize',40);
xlabel('Time (s)','FontSize',40);
subplot(212)
plot(t, vel_filtered_2)
h= gca;
set(h,'YMinorTick','off','XMinorTick','on','LineWidth',2,'Box','on','FontSize',40);
ylabel('v (m/s)','FontSize',40);
xlabel('Time (s)','FontSize',40);

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)
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
##### 1 CommentShowHide None
Quang Duy Nguyen on 19 May 2021
Hi Mathieu. Thank you so much for spending time on examining and writing the code. I've run your code and see very similar results with mine. Although we used different filters but the difference is insignificant. So, I think either the approach could be used. Cheers.