stft1.m

离散傅立叶变换用于分析吸毒病人和正常人的频谱。

fp1=fopen('e:\zyh\z07aa.dat','r');x23=fread(fp1,inf,'float');status1=fclose(fp1);
fp1=fopen('e:\zyh\b04aa.dat','r');x18=fread(fp1,inf,'float');status2=fclose(fp1);
a=20*[0 18 504 6357 47664 235494 802740 1921536 3223812 3710488 2788596 1231547 242220]-[24 684 8814 67842 347052 1241868 3184492 5891424 7799032 7200744 4399750 1596990 260388];
lp=roots(a);
b1=(exp(-(11+12*lp+3*lp.^2)./(6+11*lp+6*lp.^2+lp.^3)));
b=b1(find(b1<1&b1>0));b=b(1);
k=0:140;
w=k*2*pi;
       S1=zeros(140,141);
       SS1=zeros(140,141);
       S1(1,:)=x23(1);
       S1(2,:)=x23(2)+4*b*exp(j*w*0.008).*S1(1,:);
       S1(3,:)=x23(3)+4*b*exp(j*w*0.008).*S1(2,:)-6*b^2*exp(j*2*w*0.008).*S1(1,:);
       S1(4,:)=x23(4)+4*b*exp(j*w*0.008).*S1(3,:)-6*b^2*exp(j*2*w*0.008).*S1(2,:)+4*b^3*exp(j*3*w*0.008).*S1(1,:);
       S1(5,:)=x23(5)+4*b*exp(j*w*0.008).*S1(4,:)-6*b^2*exp(j*2*w*0.008).*S1(3,:)+4*b^3*exp(j*3*w*0.008).*S1(2,:)-b^4*exp(j*4*w*0.008).*S1(1,:);
       for n=6:140;
             S1(n,:)=x23(n)+4*b*exp(j*w*0.008).*S1(n-1,:)-6*b^2*exp(j*2*w*0.008).*S1(n-2,:)+4*b^3*exp(j*3*w*0.008).*S1(n-3,:)-b^4*exp(j*4*w*0.008).*S1(n-4,:);%公式2.5
       end
       for n=21:140;
             SS1(n-20,:)=exp(-j*w*20*0.008).*S1(n,:);
             vmin1=min(abs(SS1(n-20,:)));
             vmax1=max(abs(SS1(n-20,:)));
             SS1(n-20,:)=(abs(SS1(n-20,:))-vmin1)/(vmax1-vmin1);
       end 
        
       S2=zeros(140,141);
       SS2=zeros(140,141);
       S2(1,:)=x18(1);
       S2(2,:)=x18(2)+4*b*exp(j*w*0.008).*S2(1,:);
       S2(3,:)=x18(3)+4*b*exp(j*w*0.008).*S2(2,:)-6*b^2*exp(j*2*w*0.008).*S2(1,:);
       S2(4,:)=x18(4)+4*b*exp(j*w*0.008).*S2(3,:)-6*b^2*exp(j*2*w*0.008).*S2(2,:)+4*b^3*exp(j*3*w*0.008).*S2(1,:);
       S2(5,:)=x18(5)+4*b*exp(j*w*0.008).*S2(4,:)-6*b^2*exp(j*2*w*0.008).*S2(3,:)+4*b^3*exp(j*3*w*0.008).*S2(2,:)-b^4*exp(j*4*w*0.008).*S2(1,:);
       for n=6:140;
             S2(n,:)=x18(n)+4*b*exp(j*w*0.008).*S2(n-1,:)-6*b^2*exp(j*2*w*0.008).*S2(n-2,:)+4*b^3*exp(j*3*w*0.008).*S2(n-3,:)-b^4*exp(j*4*w*0.008).*S2(n-4,:);%公式2.5
       end
       for n=21:140;
             SS2(n-20,:)=exp(-j*w*20*0.008).*S2(n,:);
             vmin1=min(abs(SS2(n-20,:)));
             vmax1=max(abs(SS2(n-20,:)));
             SS2(n-20,:)=(abs(SS2(n-20,:))-vmin1)/(vmax1-vmin1);
       end