%REQUIRES trumpet.mat ON CD
clear;load 'trumpet.mat'
aa=100;f=800;fs=44100;L=length(X);
dt=1/fs;t=[0:L-1]*dt;
Y=X+cos(2*pi*f*t);
%GOAL: Recover X from Y.
%Design notch filter:
%Use formulae from the text.
aa2=aa*aa;wo=2*pi*f;
A=sqrt(4*aa2+aa2*aa2/wo/wo);
theta=atan(aa/wo/2);
h=-A*exp(-aa*t).*cos(wo*t+theta);
%Impulse response=delta+h.
%z=(delta+h)*y=y+h*y.
%Convolve signal and filter
%using Fourier transforms:
FY=fft(Y);Fh=fft(h)*dt;
Z=real(ifft(FY.*Fh));Z=Y+Z;
%Spectrum of trumpet+sinusoid.
%See Chapter 8 for details.
If=[1:4000];F=(If-1)*fs/L;
AY=2*abs(FY(If))/L;
AH=abs(1+Fh(If));
figure,subplot(211),plot(F,AY,F,AH,'r')
%Plot time waveforms:
It=[4601:4900];%display interval.
figure,subplot(211),plot(t(It),X(It))
figure,subplot(211),plot(t(It),Y(It))
figure,subplot(211),plot(t(It),Z(It))
figure,subplot(211),plot(t(If),h(If))
soundsc(Y,fs),pause,soundsc(Z,fs)
