signal system project

function [y,x] = IIRNotchFilter(soundfilename, notchFrequencies, alpha, beta) % y = IIRNotchFilter(soundfilename, zeroCoefficients, poleCoefficients, alpha) % % Implements a cascaded IIR digital notch filter. Each cascaded IIR notch % filter is a second order filter. Thus, for each notch frequency specified a % second order filter is added to the cascaded system. % % soundfilename A string input specifying the sound file to be processed. % notchFrequencies A vector specifying the analog frequencies to be notched. % alpha The magnitude of the poles % beta The magnitude of the betas (note alpha cannot equal beta) % % y The output processed audio sequence. % x The input audio sequence. % % Cascaded Digital Notch Filter Implementation % Dr. Ratliff [x,fs] = audioread(soundfilename); % Calculate digital frequencies of tones for each sample rate in % radians/sample and display them. w = 2*pi*notchFrequencies./fs; disp(sprintf('Digital Notch Frequencies in Rad/Sec for Fs = %d Hz',fs)); disp(w); %Calculate locations of each pole and zero, noting that each pole-zero %also has a matching complex conjugates. All I am doing here is using the %convolution command to multiply out the transfer function polynomials to %obtain the coefficients. for k=1:3 if(k==1) Z = conv([ 1 -beta*exp(1i*w(k))], [ 1 -beta*exp(-1i*w(k))]); P = conv([ 1 -alpha*exp(1i*w(k))], [ 1 -alpha*exp(-1i*w(k))]); else Z = conv(Z, [ 1 -beta*exp(1i*w(k))]); Z = conv(Z, [ 1 -beta*exp(-1i*w(k))]); P = conv(P, [ 1 -alpha*exp(1i*w(k))]); P = conv(P, [ 1 -alpha*exp(-1i*w(k))]); end end %Take the real part of the coefficients for use in the difference equation and display them B = real(Z); A = real(P); disp(sprintf('Cascaded Notch Filter Numerator (Zero) Coefficients for Fs = %d Hz',fs)); disp(B); disp(sprintf('Cascaded Notch Filter Denominator (Pole) Coefficients for Fs = %d Hz',fs)); disp(A); % Calculate the order of the difference equation O = length(A)-1; %Process the sound files using the difference equation y = zeros(size(x)); for n=O+1:length(x) y(n) = sum(x(n-O:n).*flipud(B')) + sum(y(n-O:n-1).*(flipud(-A(2:O+1)'))); end %Play Sound disp('Playing Original Sound Clip'); sound(x,fs); display('Press Any Key to Play Processed Sound Clip'); pause; sound(y,fs);