How to filter a signal that just some special frequencies can be passed?
27 次查看(过去 30 天)
显示 更早的评论
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?
0 个评论
回答(2 个)
Walter Roberson
2022-3-28
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.
3 个评论
Walter Roberson
2022-3-29
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);
Star Strider
2022-3-29
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.
.
4 个评论
Star Strider
2022-3-30
编辑:Star Strider
2022-3-30
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)
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')
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.
.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Digital Filter Design 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!