%REQUIRES trumpet.mat from CD
clear;load 'trumpet.mat'
aa=25;n=9;f=491;fs=44100;
L=length(X);dt=1/fs;t=[0:L-1]*dt;
Y=X+randn(1,L);
%GOAL: Recover X from Y.
%Note the higher noise level.
%Design comb filter:
N=poly(j*2*pi*f*[-n:n]);
D=poly(-aa+j*2*pi*f*[-n:n]);
%Compute impulse response.
%Excludes additive impulse.
%Then resonator=delta-comb
%resonator=delta-(delta+h)=-h.
[R P K]=residue(N,D);
h=-real(R.'*exp(P*t));%note the -.
figure,subplot(211),plot(t,h)
%Convolve signal and filter
%using Fourier transforms.
Fh=fft(h)*dt;FY=fft(Y);
Z=real(ifft(FY.*Fh));
%Spectrum of noisy trumpet signal:
%See Chapter 8 for details.
If=[1:4000];F=(If-1)*fs/L;
FhI=abs(Fh(If));FYI=2*abs(FY(If))/L;
figure,subplot(211),
plot(F,FhI,'r',F,FYI,'b')
%Plot time waveforms:
I=[4601:4900];%Display interval.
figure,subplot(211),plot(t(I),X(I))
figure,subplot(211),plot(t(I),Y(I))
figure,subplot(211),plot(t(I),Z(I))
soundsc(Y,fs),pause,soundsc(Z,fs)
