%REQUIRES twotrumpetsGA.mat from CD
clear;load 'twotrumpetsGA.mat'
aa=100;n=9;fr=435;fs=44100;
L=length(X);dt=1/fs;t=[0:L-1]*dt;
%Design comb filter:
N=poly(j*2*pi*fr*[-n:n]);
D=poly(-aa+j*2*pi*fr*[-n:n]);
%Compute impulse response:
%Excludes additive impulse.
[R P K]=residue(N,D);
h=real(R.'*exp(P*t));
Ih=[1:2500];%display interval
subplot(211),plot(t(Ih),h(Ih))
subplot(212),plot(t(1:300),h(1:300))
%Convolve signal and filter
%using Fourier transforms.
%Use y=(delta+h)*x=x+h*x.
Fh=fft(h)*dt;FX=fft(X);
Y=real(ifft(FX.*Fh));Y=X+Y;
%Spectrum of two-trumpet signal:
%See Chapter 8 for details.
If=[1:4000];F=(If-1)*fs/L;
FhI=abs(1+Fh(If));FXI=2*abs(FX(If))/L;
figure,subplot(211),
plot(F,FhI,'r',F,FXI,'b')
%Plot time waveforms:
figure,I=[4601:4900];%display interval
subplot(211),plot(t(I),X(I))
subplot(212),plot(t(I),Y(I))
soundsc(X,fs),pause,soundsc(Y,fs)
