c5_5_1.m

dsp入门与实践一书的源代码

clear all;
clf;
k = 0.3;  u = 1;  f=0.52;  a = 0.04;  b = 1;

samplef=60000;                 
dt=1/(10.^6);                  
sample=1/dt;                   
sampleN=1/(samplef*dt);        

N=5000;                        
F=zeros(N,1);

Zs=zeros(N,1);               	
%rand('state',sum(100*clock));
%Zs=2*rand(N,1)-1;
%Zs=b*Zs;
%plot(Zs);
%Zs=Zs-mean(Zs);

s=zeros(N,1);                 	 
x=zeros(N,1);  y=zeros(N,1);  k1=zeros(4,1);  k2=zeros(4,1);

x(1,1)=0.4;   y(1,1)=0.27;
signal=zeros(N,1);

for i=1:N
    s(i,1)=sin(2*pi*samplef*i/sample);
end
signal=a*s+Zs;

fo=45000;
while fo < 65000    
    for i=1:N
        F(i,1)=sin(2*pi*fo*i/sample);
    end

    h=2*pi*fo/sample;
    for i=2:N
        k1(1,1)=y(i-1,1);    
        k2(1,1)=x(i-1,1)-x(i-1,1).^3+f*F(i-1,1)-k*y(i-1,1)+signal(i-1,1);
         
        for j=1:2
            t1=y(i-1,1)+k2(j,1)*h/2;
            t11=x(i-1,1)+k1(j,1)*h/2;
            t3=(F(i-1,1)+F(i,1))/2;
            t4=(signal(i-1,1)+signal(i,1))/2;        
            k1(j+1,1)=t1;
            k2(j+1)=t11-t11.^3+f*t3-k*t1+t4;
        end
  
        t2=y(i-1,1)+k2(3,1)*h;
        t22=x(i-1,1)+k1(3,1)*h;
        k1(4,1)=t2;
        k2(4,1)=t22-t22.^3+f*F(i,1)-k*t2+signal(i,1);   
        x(i,1)=x(i-1,1)+(k1(1,1)+2*k1(2,1)+2*k1(3,1)+k1(4,1))*h/6;
        y(i,1)=y(i-1,1)+(k2(1,1)+2*k2(2,1)+2*k2(3,1)+k2(4,1))*h/6;
    end
   
    zpoint=zeros(N,1);
    m=1;

    for i=1:N-1
        if x(i,1)==0
            zpoint(m,1)=i;
            m=m+1;
        elseif x(i,1)*x(i+1,1) < 0
            zpoint(m,1)=i;
            m=m+1;
        end
    end

     m=m-2;
    zspace=zeros(m,1);
    for i=1:m
        zspace(i,1)=zpoint(i+1,1)-zpoint(i,1);
    end

    m=m-1;
    difzs=zeros(m,1);
    
    for i=1:m
        difzs(i,1)=abs(zspace(i+1,1)-zspace(i,1));
    end

    result=zeros(m,1);
    i=1;    j=1;

    while i < m
        counter=0;
        while difzs(i,1) <= 3 
            if i>m-1
                break;
            else
                i=i+1;
            end
        end
  
        while difzs(i,1)> 3
            if i>m-1
                break;
            else
                i=i+1;
            end 
        end
  
        Nu=i;
        while difzs(i,1) <= 3 
            counter=counter+1;
            if i>m-1
                break;
            else
                i=i+1;
            end 
            if counter>30
                result(j,1)=zpoint(Nu,1);
                j=j+1;
                counter=0;
                break;
            end
        end
    end
    j=j-1;
    
    dd=zeros(j,1);
    ddd=zeros(j,1);

    for i=1:j-1
        dd(i,1)=result(i+1,1)-result(i,1);
    end

    for i=1:j-2
        ddd(i)=abs(dd(i+1,1)-dd(i,1));
    end

    for i=1:j-2
        if ddd(i,1)> 250;
            for k=1:j
                result(k,1)=0;
            end
                break;
        end
    end

    if j==1 | j==0
        avg=0;
    else
        for i=1:j
            avg = result(j,1)-result(1,1);
            avg=avg/(j-1);
        end
    end
    
    if avg > 0 & fo < samplef
        avg1=avg;
        fo=fo*1.03;
    elseif avg > 0 & fo > samplef
        fo=45000;
        break;
    elseif avg ==0 & fo > samplef
        fo=45000;
        break;
    else
        fo=fo*1.03
    end
end