T=1/8000;
f0=300;
k=0:15999;
x=cos(2*pi*f0*k*T+0.2);

h=1/T*[1 -2*cos(2*pi*f0*T) 1];

y=T*conv(x,h);


f=1/T*(0:511)/512;

X=zeros(1,512);

for n=1:length(h),
    X=X+h(n)*exp(-j*2*pi*f*(n-1)*T);
end;
X=T*X;

plot(f,abs(X))


x=cos(2*pi*f0*k*T+0.2)+cos(2*pi*200*k*T+0.2)+cos(2*pi*3000*k*T+0.2);
y=T*conv(x,h);


load frase_8192.mat;
orig=orig/max(abs(orig));
k=0:length(orig)-1;

xf=orig+cos(2*pi*f0*k*T+0.2);
y=T*conv(xf,h);
soundsc(y,8000);


C=sum([1 -2*0.98*cos(2*pi*f0*T) 0.98^2])/sum([1 -2*cos(2*pi*f0*T) 1]);

b=C/T*[1 -2*cos(2*pi*f0*T) 1];
a=[1 -2*0.98*cos(2*pi*f0*T) 0.98^2];

[H,F] = freqz(b,a,512,8000);
plot(F,T*abs(H));
y=filter(b,a,xf);
soundsc(y,8000);