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

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)

Answer 1

[Star Strider]

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.

[Manuel Cruz]

Thank you, just one last question, where do you take the Rp and Rs values? To add it to my bibliography.

[Star Strider]

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.

[Manuel Cruz]

I really appreciate your help, thanks.

[Star Strider]

As always, my pleasure!