%-----------------------------------
% script spectrum
clear
%load handel
%sound(y,Fs)
[y,Fs,bits]=wavread('tele_1.wav');
Fs
a=input('starting point')
n=input('length of series (=power of 2 is fastest)')
x=y(a:a+n); m=ceil(n/2); 
sound(x,Fs)
%-----------------------------------
z = fft(x(1:n));
% create a structure P2{a,b}
% P2.a is a vector of size m+1 (cosine coefficients)
% P2.b is a vector of size m-1 (sine coefficients)
% to find the amplitude associated with frequency n
% we must sum (in the pythagorean sense) its cosine and
% sine components
P2 = struct('a',real(z(1:m+1))/m,'b',-imag(z(2:m))/m);
P2.a(1)=P2.a(1)/2;P2.a(m+1)=P2.a(m+1)/2;
E(2:m)=sqrt(P2.a(2:m).^2+P2.b(1:m-1).^2);
E(1)=abs(P2.a(1));E(m+1)=abs(P2.a(m+1));
Bvals = linspace(0,n/Fs,n+1)';
%Gvals = interpft(z);
%Gvals = CSeval(P2,n/Fs,Bvals);
subplot(2,1,1);
plot(Bvals,x)
subplot(2,1,2);
plot(E(1:200))
for i=1:m+1
    if (E(i)>.038)
        i-1
    end
end
%axis([0 360 -.05 1.05])
%set(gca,'xTick',linspace(0,360,13))