% Demo Transformada Chirp
% Primero leemos la señal correspondiente a do4.wav. Es
% un segmento de 256 muestras de do4 ... 261.6 Hz.

echo on;
[x,Fm,Nbits]=wavread('do4.wav');
figure(1);
clf;
plot(x);
title('forma de onda, do 4, frec. fundamental= 261.6 Hz');
xlabel('muestras');
ylabel('amplitud');
grid on;
pause

% calculamos su TF con 1024 puntos

tic
X=fft(x,2^10);
toc
clf;
f0=0;
ff=(2^9-1)*Fm/2^10;
df=Fm/2^10;
f=[f0:df:ff];
plot(f,abs(X(1:2^9)));
hold on
plot(261.6,1.78,'m*');
title('DFT de 2\^10 puntos');
xlabel('frecuencia analógica');
ylabel('amplitud');
grid on;
pause
clf
stem(f,abs(X(1:2^9)));
hold on
plot(261.6,1.78,'m*')
axis([230 290 0 1.85]);
xlabel('frecuencia analógica');
ylabel('amplitud');
grid on
title('Detalle de la DFT alrededor de la frec. fundamental');
pause

% calculamos un zoom sobre la frecuencia fundamental de 261.6 Hz

f0=240;
W=exp(-j*2*pi/2^18);
A=exp(j*2*pi*f0/Fm);
M=1024;
tic
for i=1:10
   G=czt(x,M,W,A);
end   
toc/10
hold on
df=1/2^18*Fm;
ff=f0+df*(M-1);
plot([f0:df:ff],abs(G),'r');
axis([230 290 0 1.85]);
[ymax,xmax]=max(abs(G));
fmax=f0+df*xmax;
plot(fmax,ymax,'r*');
title(['En rojo el resultado de 1024 puntos Transf. Chirp. Frec. max ' num2str(fmax)]);

pause

% calculamos el zoom con la FFT de 2^18 puntos= 262144 puntos

tic
for i=1:10
   XF=fft(x,2^18);
end   
toc/10
f0=0;
ff=(2^17-1)*Fm/2^18;
df=Fm/2^18;
f=[f0:df:ff];
plot(f,abs(XF(1:2^17)),'g');
[ymax,xmax]=max(abs(XF(1:6900))); %buscamos sobre le primer armónico
fmax=df*xmax;
plot(fmax,ymax,'g*');
title(['En verde DFT de 2\^18 (262144) puntos. Frec. max ' num2str(fmax)]);
echo off