goertzel算法.m

GOERTZEL算法,应用于音调检测

% MATLAB 设计GOERTZEL算法
clear all;
tab=[1.7077;1.6453;1.5687;1.4782;1.1641;0.9964;0.7986;0.5685];
tab4=[-0.52047;-0.56857;-0.62033;-0.67358;-0.81314;-0.86706;-0.91680;-0.95874];
number=input('number=','s');
w=tab;
k=0;
m=0;
if (number=='0')
    k=4;
    m=6;
elseif(number=='1')
    k=1;
    m=5;
elseif(number=='2')
    k=1;
    m=6;
elseif(number=='3')
    k=1;
    m=7;
elseif(number=='4')
    k=2;
    m=5;
elseif(number=='5')
    k=2;
    m=6;
elseif(number=='6')
    k=2;
    m=7;
elseif(number=='7')
    k=3;
    m=5;
elseif(number=='8')
    k=3;
    m=6;
elseif(number=='9')
    k=3;
    m=7;
elseif(number=='a')
    k=1;
    m=8;
elseif(number=='b')
    k=2;
    m=8;
elseif(number=='c')
    k=3;
    m=8;
elseif(number=='d')
    k=4;
    m=8;
elseif(number=='*')
    k=4;
    m=5;
elseif(number=='#')
    k=4;
    m=7;
end
    y1(1)=tab4(k);
    y1(2)=0;
    y2(1)=tab4(m);
    y2(2)=0;
for i=3:205
    y1(i)=tab(k)*y1(i-1)-y1(i-2);
    y2(i)=tab(m)*y2(i-1)-y2(i-2);
    x(i)=y1(i)+y2(i);
end
%bx=max(abs(x))
%zq=quantizer('fixed', 'ceil', 'saturate',[16 15]);
%a=num2hex(zq,y1)
%a=y1*2^15;
figure(1);
subplot(311);plot(x);
xx=fft(x,205);
pyy=xx.*conj(xx);
f=(0:(205/2-1));
for i=1:205/2-1
    f(i)=f(i)*8000/205;
end
subplot(312);plot(f,pyy(1:205/2));
tab2=[1.7033;1.6359;1.5623;1.4829;1.1631;1.0088;0.7901;0.5595]; %建立表矩阵1次谐波 N=205 coeff 
tab3=[0.9011;0.6760;0.4408;0.1989;-0.6471;-0.9823;-1.3758;-1.6870];
vk(1)=0;
vk(2)=0;
w=tab2;
for i1=1:8
    for i2=3:205
        vk(i2)=w(i1)*vk(i2-1)-vk(i2-2)+x(i2);
    end
    a=vk(205)^2+vk(204)^2;
    b(i)=w(i1)*vk(205);
    c(i)=b(i)*vk(204);
    Xk(i1)=vk(205)^2+vk(204)^2-w(i1)*vk(205)*vk(204);
end

for i=1:3
    if(Xk(i+1)<Xk(i))
        if(Xk(i)>7000)
            sq(1)=i;
        end
    elseif(Xk(i+1)>7000)
        sq(1)=i+1;
    end
end
for i=5:7
    if(Xk(i+1)<Xk(i))
        if(Xk(i)>7000)
            sq(2)=i;
        end
    elseif(Xk(i+1)>Xk(i)&&Xk(i+1)>7000)
        sq(2)=i+1;
    end
end
w2=tab3;
    for i2=3:205
        vk(i2)=w2(sq(1))*vk(i2-1)-vk(i2-2)+x(i2);
    end
    Xk2(1)=vk(205)^2+vk(204)^2-w2(sq(1))*vk(205)*vk(204);

    for i2=3:205
        vk(i2)=w2(sq(2))*vk(i2-1)-vk(i2-2)+x(i2);
    end
    Xk2(2)=vk(205)^2+vk(204)^2-w2(sq(2))*vk(205)*vk(204);
str='123A456B789C*0#D';
q=(sq(1)-1)*4+sq(2)-4;
re=str(q);
%figure(2);
subplot(313);stem(Xk);title(['你输入的号码是',re],'Color','r');