MATLAB Answers

Problems using a Butterworth filter instead of a bandstop to filter frequencies.

1 view (last 30 days)
Manuel Cruz
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

Sign in to comment.

Accepted Answer

Star Strider
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

Show 1 older comment
Star Strider
Star Strider on 29 Nov 2020
As always, my pleasure!
I just use the values for them that are likely appropriate. They can be anything reasonable, however anything less tha 1 dB (passband attenutation or ripple amplitude) or greater than 100 dB (stopband attenuation or ripple amplitude) can be next to impossible to realise in any filter, so I choose something in between that would appear to give reasonable results in a particular context. Everything depends on the signal and the application.
These are design parameters, so there is no straightforward way to estimate them from the data (however if you had input and output vectors and you wanted to estimate the transfer function that produced them, the System Identification Toolbox would likely provide reasonable estimates.)
Obviously in this application, use the value for ‘Fs’ from the file. So long as it is greater than about 365 Hz, the filter will work.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!