# How to filter a signal that just some special frequencies can be passed?

40 views (last 30 days)
Mitra Taghizadeh on 28 Mar 2022
Edited: Star Strider on 30 Mar 2022
Hello every one.
I have a random vibration signal. I want to filter it in the way that just three special frequencies passed and the other recorded data should be filtered. Do you know how can I do that?

Walter Roberson on 28 Mar 2022
fft() it with the number of points set such that all three frequencies correspond to integer multiples of the sampling frequency. Then zero out all of the bins except those ones (and their complex conjugects), then ifft() back
Using traditional filter banks will not work well: they need to ramp between "present" and "not-present" so they extract a smoothed version of the frequency components instead of just those frequency components.
Mitra Taghizadeh on 29 Mar 2022
I have done the work that you say. but there is another problem.
I have 6400 data. After applying FFT it decreases to 3200. After filtering prosidure and applying IFFt, I will have 3200 data which is not comparable to the unfiltered data because of size. When I put them in a figure it seems that there is some phase difference.

Star Strider on 29 Mar 2022
A reasonably efficient way of doing this is to use a multiband FIR filter, for example —
Fs = 1000; % Use Correct Sampling Frequency (Must Be Greater Than 370 Hz)
fcomb = [[55 59 61 64]+60, [55 59 61 64]+180, [55 59 61 64]+240]; % Band Frequencies (Hz)
mags = [[0 1 0] [1 0] [1 0]]; % Magnitudes Of Passbands
% mags = [[1 0 1], [0 1], [0 1]]; % Magnitudes Of Stopbands
dev = [[0.1 0.5 0.1], [0.5 0.1], [0.5 0.1]]; % Passband Deviation Allowances
% dev = [[0.5 0.1 0.5], [0.1 0.5], [0.1 0.5]]; % Stopband Deviation Allowances
[n,Wn,beta,ftype] = kaiserord(fcomb,mags,dev,Fs);
hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');
figure
freqz(hh, 1, 2^20, Fs)
set(subplot(2,1,1), 'XLim', [0 500]) % Zoom X-Axis
set(subplot(2,1,2), 'XLim', [0 500]) % Zoom X-Axis % DataFilt = filtfilt(hh, 1, Data); % Filter Signal
The disadvantage of a multiband FIR filter such as this is that it can require a very long signal, although using it with the fftfilt function may offer a way around that. If that is not an option, a series of parallel IIR filters (specifically elliptic filters, since they are the most computationally efficient) may be necessary. The approach depends on the signal itself and the signal processing constraints.
.
Star Strider on 30 Mar 2022
Edited: Star Strider on 30 Mar 2022
The FIR filter is too long, nearly 1.7 times the length of the signal. Since it has to be less than half the signal length, it will be necessary to use IIR filters in parallel —
Fs = 128; % Use Correct Sampling Frequency (Must Be Greater Than 370 Hz)
% fcomb = [[55 59 61 64]+60, [55 59 61 64]+180, [55 59 61 64]+240]; % Band Frequencies (Hz)
fcomb = [[0.37 0.39 0.43 0.44] [0.44 0.45 0.47 0.49] [0.58 0.61 0.63 0.65] [1.58 1.60 1.64 1.68]];
mags = [[0 1 0] [1 0] [1 0] [1 0]]; % Magnitudes Of Passbands
% mags = [[1 0 1], [0 1], [0 1]]; % Magnitudes Of Stopbands
dev = [[0.1 0.5 0.1], [0.5 0.1], [0.5 0.1] [0.5 0.1]]; % Passband Deviation Allowances
% dev = [[0.5 0.1 0.5], [0.1 0.5], [0.1 0.5]]; % Stopband Deviation Allowances
[n,Wn,beta,ftype] = kaiserord(fcomb,mags,dev,Fs);
hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');
FilterLength = numel(hh)
FilterLength = 10745
figure
freqz(hh, 1, 2^20, Fs)
set(subplot(2,1,1), 'XLim', [0 2]) % Zoom X-Axis
set(subplot(2,1,2), 'XLim', [0 2]) % Zoom X-Axis Fs = 128; % Sampling Frequency
Fn = Fs/2; % Nyquist Frequency
Rs = 50; % Stopband Ripple/Attenuation
Rp = 1; % Passband Ripple/Attenuation
Wp = [0.39 0.43; 0.44 0.48; 0.60 0.64; 1.60 1.64]; % Passband Matrix (Hz)
for k1 = 1:size(Wp,1)
Wpn = Wp(k1,:)/Fn; % Normalised Passband For This Filter
Ws = Wpn.*[0.95 1.05]; % Normalised Stopband For This Filter
[n,Wn] = ellipord(Wpn,Ws,Rp,Rs); % Order Estimation
[z,p,k] = ellip(n,Rp,Rs,Wn); % Zero-Pole-Gain Filter Realisation
[sos{k1},g(k1)] = zp2sos(z,p,k); % Second-Order-Section Implementation For Stability
figure
freqz(sos{k1}, 2^18, Fs)
set(subplot(2,1,1), 'XLim',[0 2]) % Zoom X-Axis
set(subplot(2,1,2), 'XLim',[0 2]) % Zoom X-Axis
sgtitle(sprintf('Bandpass: %.2f Hz to %.2f Hz',Wp(k1,:)))
end    M1 = 6400×2
0.0078 0.6030 0.0156 -0.6030 0.0234 0.4020 0.0312 0.6030 0.0391 0.4020 0.0469 0 0.0547 0 0.0625 -0.4020 0.0703 0.4020 0.0781 0.4020
DataFilt = zeros(size(M1,1), size(Wp,1)); % Preallocate
for k1 = 1:size(Wp,1)
DataFilt(:,k1) = filtfilt(sos{k1}, g(k1), M1(:,2)); % Filter Signal
end
figure
subplot(2,1,1)
plot(M1(:,1), M1(:,2))
grid
title('Unfiltered Signal')
subplot(2,1,2)
plot(M1(:,1), sum(DataFilt,2))
grid
title('Filtered Signal')
sgtitle('Results') y = M1(:,2);
FTy = fft(y-mean(y))/numel(y); % Calculate & Plot Fourier Transform
Fv = linspace(0, 1, numel(y)/2+1)*Fn;
Iv = 1:numel(Fv);
figure
plot(Fv, abs(FTy(Iv))*2)
grid
xlim([0 2])
xlabel('Frequency')
ylabel('Amplitude')
title('Original Signal Spectrum') y = sum(DataFilt,2);
FTy = fft(y-mean(y))/numel(y); % Calculate & Plot Fourier Transform
Fv = linspace(0, 1, numel(y)/2+1)*Fn;
Iv = 1:numel(Fv);
figure
plot(Fv, abs(FTy(Iv))*2)
grid
xlim([0 2])
xlabel('Frequency')
ylabel('Amplitude')
title('Filtered Signal Spectrum') The elliptic IIR filters and the parallel implementation of them appear to have performed as desired!
EDIT — (30 Mar 2022 at 11:10)
Expanded comment documentation. Code unchanged.
.