Notch Filter Using Z-Transform Equations and Filter Functions Only

I am provided a signal with two frequency disturbance peaks and I need to develop a notch filter using the 'filter' function in MATLAB to remove each. I know the frequencies I want to remove (1200 Hz and 2500 Hz), and so I developed the Z-transform functions

(Z - 0.999*exp(1j*0.054422*pi))(Z - 0.999*exp(-1j*0.054422*pi))/...

(Z - 0.99*exp(1j*0.054422*pi))(Z - 0.99*exp(-1j*0.054422*pi))

and

(Z - 0.999*exp(1j*0.113379*pi))(Z - 0.999*exp(-1j*0.113379*pi))/...

(Z - 0.99*exp(1j*0.113379*pi))(Z - 0.99*exp(-1j*0.113379*pi))

I then found the coefficients of difference equations and included them into the code as arrays 'b', 'a', 'b1' and 'a1'.

My issue is that my output of filt_2 looks no different than the original signal. Overall my filters looks exactly as I would expect. My code and images below.

Please advise on why my filter portion is not filtering.

%Read in audio and convert to frequency domain
[y Fs] = audioread('my_audio_1.wav');
N = length(y);
n = 0:N-1;
OM = 0:0.01:pi;
X = exp(-1j*OM'*n)*y;
fq = Fs*OM/(2*pi);

%filt_1 coefficients for removing 1200 Hz disturbance
b = [1 -0.999*(exp(-1j*0.054422*pi)+exp(1j*0.054422*pi)) 0.998001];
a = [1 -0.99*(exp(-1j*0.054422*pi)+exp(1j*0.054422*pi)) 0.9801];
%filt_2 coefficients for removing 2500 Hz disturbance
b1 = [1 -0.999*(exp(-1j*0.113379*pi)+exp(1j*0.113379*pi)) 0.998001];
a1 = [1 -0.99*(exp(-1j*0.113379*pi)+exp(1j*0.113379*pi)) 0.9801];

%cascaded filters
filt_1 = filter(b,a,X); %filter the original signal
filt_2 = filter(b1,a1,filt_1); %filter the first filtered signal

%output
tiledlayout(2,1);
nexttile;
plot(fq, abs(X));
nexttile;
plot(fq, abs(filt_2),'r');

figure();
zplane(b,a);
hold;
zplane(b1,a1);
fvtool(b,a);
fvtool(b1,a1);

Answer 1

[Star Strider]

Your filter works correctly, although your function for analysing and plotting the filtered result fails.

Try this:

%filt_1 coefficients for removing 1200 Hz disturbance
b = [1 -0.999*(exp(-1j*0.054422*pi)+exp(1j*0.054422*pi)) 0.998001];
a = [1 -0.99*(exp(-1j*0.054422*pi)+exp(1j*0.054422*pi)) 0.9801];
%filt_2 coefficients for removing 2500 Hz disturbance
b1 = [1 -0.999*(exp(-1j*0.113379*pi)+exp(1j*0.113379*pi)) 0.998001];
a1 = [1 -0.99*(exp(-1j*0.113379*pi)+exp(1j*0.113379*pi)) 0.9801];

Fn = 2.5E+4;                    % Nyquist Frequency (Hz)
Fs = Fn*2;                      % Sampling Frequency (Hz)

[h,f] = freqz(b,a,2^14,Fs);
[h1,f1] = freqz(b1,a1,2^14,Fs);

figure
subplot(2,1,1)
plot(f, 20*log10(abs(h)))
hold on
plot(f1, 20*log10(abs(h1)))
hold off
xlabel('Frequency (Hz)')
ylabel('Amplitude (dB)')
grid
subplot(2,1,2)
plot(f, rad2deg(unwrap(angle(h))))
hold on
plot(f1, rad2deg(unwrap(angle(h1))))
hold off
xlabel('Frequency (Hz)')
ylabel('Phase (°)')
grid

y = randn(1,Fs)';
N = length(y);

%cascaded filters
filt_1 = filter(b,a,y); %filter the original signal
filt_2 = filter(b1,a1,filt_1); %filter the first filtered signal

FTy = fft(y)/N;
FTfilt_2 = fft(filt_2)/N;
Fv = linspace(0, 1, fix(N/2)+1);
Iv = 1:numel(Fv);

TrFn = FTfilt_2 ./ FTy;         % Filter Transfer Function

%output
figure
tiledlayout(2,1);
nexttile;
plot(Fv, abs(FTy(Iv)));
nexttile;
plot(Fv, abs(TrFn(Iv)),'r');

I have no idea what you intended with ‘X’, however it is not performing any sort of Fourier transform that I recognise.