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

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?

Answer 1

[Walter Roberson]

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]

Thank you for your response.

How can I zero out all of the bins? Could you please help me?

[Walter Roberson]

to_save = [17 84 90]  %bin numbers of frequencies to save
temp = zeros(size(FFT_result));
temp(to_save) = 1;
temp(end - to_save + 2) = 1;
selected_FFT = FFT_result .* temp;
filtered_signal = ifft(selected_FFT);

[Mitra Taghizadeh]

Thanks for your help.

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.

Would you please help me for this issue?

Answer 2

[Star Strider]

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.

.

[Mitra Taghizadeh]

Thank you for your response.

Can I ask you to give a bief explenation about that how we should choose "mags". You know in my code the perposed frequencies are close to each other and near zero (these frequecies should pass: .14 .46 1.6 ). could you help me to define appropriate fcomb mags and dev. Becuse I got some errors from MATLAB

[Star Strider]

My pleasure.

See the documentation on kaiserord for a full discussion.

The ‘mags’ (called a in the documentation) vector values are either 0 or 1 here although they can apparently be anything. They define the amplitudes in the vector (‘comb’ here) that define the band edges. So in the first filter segment, the frequency band defined by [55 59] has 0 amplitude, the band defined by [59 61] has an amplitude of 1 and so for the rest.

The ability to even design the filter depends on what the sampling frequency is (that is still unknown), and how long the time-domain vector is that you want to filter. (For frequencies as close together as stated, a FIR filter may not be appropriate. Three IIR elliptic filters in parallel would likely be required, and depending on the sampling frequency, even those might be a challenge to design. I will not be able to determine that without knowing ther vector length and the sampling frequency.)

[Mitra Taghizadeh]

Thanks for your reply.

I have attached the data. Sampling frequency is 128 Hz. and I want to pass just dominant frequencies including (0.41 0.46 0.62 1.62).

[Star Strider]

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 = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/945529/data.xlsx')

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.

.