Problems using a Butterworth filter instead of a bandstop to filter frequencies.
29 views (last 30 days)
Show older comments
Manuel Cruz
on 29 Nov 2020
Commented: Star Strider
on 29 Nov 2020
I'm building a bandstop filter but the bandstop function takes a long long time to run, so I decided to try filtering with a Butterworth filter, but I don't get the same response as with bandstop. Can you help me please?
%% Frequency range
lfreq = 170; % Lower freq
hfreq = 180; % Highest freq
%% Extract audio signal
[f Fs] = audioread('Schumann.wav'); % Extract signal
info = audioinfo('Schumann.wav'); % Extract audio info
n = length(f); % Length of signal
t = 0:seconds(1/Fs):seconds(info.Duration); % Create a x-axis of time in s
t = t(1:end-1); % Set time
freq = 1/time2num(t(2)-t(1))/n*(0:n); % Create a x-axis of frequencies in Hz
freq = freq(1:end-1); % Set frequency
L = 1:floor(n/2); % Only take the first half of freq
%% Compute the Fast Fourier Transform
fhat = fft(f,n); % Compute the FFT
PSD = fhat.*conj(fhat)/n; % Power spectum (energy per freq)
%% Use "bandstop" to filter frequency range
% ffilt = bandstop(f,[lfreq hfreq],Fs); % Remove all freq in range
% fhatfilt = fft(ffilt,n); % Compute the FFT for filtered signal
% PSDfilt = fhatfilt.*conj(fhatfilt)/n; % Filtered power spectum
%% Butterworth filter <----
Wn = [lfreq hfreq]/n;
[b a] = butter(5,Wn,'stop');
ffilt = filtfilt(b,a,f);
%% Plots (original and filtered signal)
figure(1)
plot(freq(L),PSD(L),'r','LineWidth',1.5); hold on
plot(freq(L),PSDfilt(L),'b','LineWidth',1.5);
legend('Original','Filtered');
xlabel('Frequency [Hz]'); ylabel('PSD [J/Hz]');
figure(2)
plot(time2num(t),f,'r','LineWidth',1.5); hold on
plot(time2num(t),ffilt,'b','LineWidth',1.5);
legend('Original','Filtered');
xlabel('Time [s]'); ylabel('Amplitude [Unknown]');
%sound(ffilt,Fs)
0 Comments
Accepted Answer
Star Strider
on 29 Nov 2020
You are not designing the filter correctly.
Try this instead:
lfreq = 170; % Lower freq
hfreq = 180; % Highest freq
Fs = 44100; % Sampling Frequency
Fn = Fs/2; % Nyquist Frequency
Ws = [lfreq hfreq]/Fn; % Stopband Frequency (Normalised)
Wp = [0.99 1.01].*Ws; % Passband Frequency (Normalised)
Rp = 1; % Passband Ripple
Rs = 90; % Passband Ripple (Attenuation)
[n,Wp] = ellipord(Wp,Ws,Rp,Rs); % Elliptic Order Calculation
[z,p,k] = ellip(n,Rp,Rs,Wp,'stop'); % Elliptic Filter Design: Zero-Pole-Gain
[sos,g] = zp2sos(z,p,k); % Second-Order Section For Stability
figure
freqz(sos, 2^18, Fs) % Filter Bode Plot
set(subplot(2,1,1), 'XLim',Wp*Fn) % Optional
set(subplot(2,1,2), 'XLim',Wp*Fn) % Optional
Use the filtfilt function to do the actual filtering.
4 Comments
More Answers (0)
See Also
Categories
Find more on Digital Filter Analysis in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!