%REQUIRES twotrumpetsAB.mat from CD
clear;load 'twotrumpetsAB.mat'
fr=436;%Note A is to be rejected
fs=44100;aa=0.99;N=length(X);
n=round(fs/fr);%See Comments
B=[1 zeros(1,n) -1];
A=[1 zeros(1,n) -aa^n];
W=2*pi/N*[0:4999];EW=exp(j*W);
H=polyval(B,EW)./polyval(A,EW);
FX=10*abs(fft(X))/N; subplot(211),
plot(W,FX(1:5000),W,abs(H),'r')
Y=filter(B,A,X);
I=[1000:1199];T=I/fs;figure
subplot(211),plot(T,X(I))
subplot(212),plot(T,Y(I))
soundsc(X,fs),pause,soundsc(Y,fs)
