演示程序.zip

clear all; close all; clc; tic; x=[10 01 11 00 10 11 01 00 10 01]; b=5; N=length(x); Ts=0.00001; fs=1/Ts; %%100KHZ Tm=200*Ts; EndTime=2*Tm-Ts; %跳频周期0.004s Nset=(EndTime+Ts)/Ts;%%400 y=zeros(Nset,1); i=0; xd=zeros(1,N); HTT=[];HTT1=[];HTT2=[]; while i<N i=i+1; switch (x(i)) case 00 [u1,time1]=gensig('sin',Tm/1,EndTime,Ts); xd(i)=1; case 01 [u1,time1]=gensig('sin',Tm/2,EndTime,Ts); xd(i)=2; case 10 [u1,time1]=gensig('sin',Tm/3,EndTime,Ts); xd(i)=3; case 11 [u1,time1]=gensig('sin',Tm/4,EndTime,Ts); xd(i)=4; otherwise disp( '输入错误'); end y(1+(i-1)*Nset:Nset*i,1)=u1; %t=[t;time1]; end %跳频信号时域波形 figure(1) subplot(2,1,1); plot(y) xlabel('t') ylabel('4FSK调制信号') title('4FSK信号') %axis([0 20 -1 1]) grid %输出跳频序列1111111111111111111111111111111 G=[0 1 2 3 4 5 6 7 8 9 10]; Su=mod(b.*G,11); for j=2:11 fprintf('%3d',Su(j)) end Tc=30*Ts;%%%3KHZ--30KHZ,10个频点，带宽3KHZ ff1=[]; C0=zeros(Nset,1); for k=2:(N+1) TT=1/(Su(k)/Tc+xd(k-1)/Tm); [Carrier1,time2]=gensig('sin', TT,EndTime,Ts); % ff1=Su(k)/Tc; % ff=[ff ff1];%%载频 C0(1+(k-2)*Nset:Nset*(k-1),1)=Carrier1; end %混频产生跳频信号 subplot(2,1,2) plot(C0) xlabel('t') ylabel('跳频调制信号') title('跳频信号') %%%%%%%%%%%%%%修正%%%%%%%%%%%%%%%%%%%%%% NT=EndTime*fs+1; cc1=0; ae=1; ER1=zeros(1,1); ER2=zeros(1,1); for SNR=20:20 j1=0; errsum1=0; ersum=0; ave_peroid=[];ave_peroid1=[];ave_peroid2=[]; N1=1; while(j1<N1) FHsignal=awgn(C0,SNR);%加噪 %%%%%%%%%%%时频分析%%%%%%%%%%%%%% z=hilbert(FHsignal); e=2^nextpow2(length(z)/4); % figure(2); h1= WINDOW(401,'Hamming'); [tfr,t,f]= TFRSTFT(z,1:length(z),e,h1,0); %% d=length(f); t1=length(t); cc=f(1:d/2); fmax0=max(abs(tfr)); fFH0=abs(fft(fmax0)); [maxr0,f0]=max(abs(fFH0(2:length(fFH0)))); NN0=length(fmax0)/f0; er=abs(NN0-NT)/NT; ersum=ersum+er; %STFT的误差 figure(3) contour(t,cc*2,abs(tfr)); axis([0 4000 0 0.5]) A=abs(tfr).^2; figure(4) contour(t,cc*2,A); axis([0 4000 0 0.5]) B=[]; Amax=max(max(A)); for q=1:t1 A1=A(:,q); %取列向量 A3=[];Z=[];D2=[];Y1=0;C=64;L=16; for c=1:C A2=A1((c-1)*L+1:c*L); A3=[A3 A2]; % 分段后矩阵元素的集合 Y(c)=sum(A2); end r0=Amax*0.6; A4=[]; for c=1:C if Y(c)>r0 %信号段 A3(:,c)=1; else %噪声段 A3(:,c)=0; end end for c=1:C A4((c-1)*L+1:c*L)=A3(:,c)'; end B=[B A4']; %得到二值化后的时频矩阵 end % 匹配搜索算法 BP=[]; for p=1:d/2 bp=B(p,:); %取B的行向量 ap=A(p,:); %取A的行向量 bm=[];M=0; for r=1:q-1 if bp(r)~=bp(r+1) bm=[bm r+1]; M=M+1;%把发生变化的起点标注出来 end end bpm=[]; if M==0 %全0的情况 else fm=[]; for m=1:M-1 bpm=bp; for n=bm(m):bm(m+1)-1 bpm(n)=mod(bp(n)+1,2);%对第m组取反 end fbp_hdmd=[];fbp1=[];fbpD=[0];k1=0;k2=0;fbp=0; fbpm_hdmd=[];fbpm1=[];fbpmD=[0];km1=0;km2=0;fbpm=0; uu1=0.5;uu2=0.6; %调这两个参（uu1在0.4~0.6，uu2在0.4~0.8） for r1=1:q fbp_hdmd=[fbp_hdmd bp(r1)*ap(r1)]; fbpm_hdmd=[fbpm_hdmd bpm(r1)*ap(r1)]; end for r2=1:q fbp1(r2)=ap(r2)-fbp_hdmd(r2); fbpm1(r2)=ap(r2)-fbpm_hdmd(r2); end for r3=2:q fbpD(r3)=bp(r3)-bp(r3-1); fbpmD(r3)=bpm(r3)-bpm(r3-1); end for r4=1:q if(bp(r4)~=0) k1=k1+1; end if(fbpD(r4)~=0) k2=k2+1; end if(bpm(r4)~=0) km1=km1+1; end if(fbpmD(r4)~=0) km2=km2+1; end end fbp=sum(fbp1)+uu1*Amax*k1+uu2*Amax*k2; fbpm=sum(fbpm1)+uu1*Amax*km1+uu2*Amax*km2; fm=[fm fbpm-fbp]; end if(min(fm)<0) [~,hh]=find(min(fm)); for n1=bm(hh):bm(hh+1)-1 bpm(n1)=mod(bp(n1)+1,2);%对第m组取反 end bp=bpm; end end BP=[BP bp']; end BPP=BP'; figure(5) contour(t,cc*2,BPP); axis([0 4000 0 0.5]) %传统算法 [fmax,ff]=max(A); fFH1=abs(fft(fmax)); [maxr1,f1]=max(fFH1(2:length(fFH1))); EN1=length(z)/f1;%长度估计 %subplot(3,1,3) ht=ff(2:length(ff)); ht1=ff(1:length(ff)-1)-ht; htt2=abs(ht1(1:length(ht1))); htt2=htt2/norm(htt2); figure(10) contour(t,cc*2,A); hold on plot(htt2) axis([0 4000 0 0.5]) hold off a1=sort(htt2,'descend'); vv1=find(htt2>0.3*max(htt2)); peroid1=[]; for l=1:length(vv1)-1 peroid1=[peroid1 vv1(l+1)-vv1(l)]; end ave_peroid1=[ave_peroid1 sum(peroid1)/(length(vv1)-1)]; %优化算法 [fmax,ff]=max(BPP); fFH1=abs(fft(fmax)); [maxr1,f1]=max(fFH1(2:length(fFH1))); EN1=length(z)/f1;%长度估计 %subplot(3,1,3) ht=ff(2:length(ff)); ht1=ff(1:length(ff)-1)-ht; htt2=abs(ht1(1:length(ht1))); htt2=htt2/norm(htt2); figure(11) contour(t,cc*2,BPP); hold on plot(htt2) axis([0 4000 0 0.5]) hold off a=sort(htt2,'descend'); vv=find(htt2>0.3*max(htt2)); peroid=[]; for l=1:length(vv)-1 peroid=[peroid vv(l+1)-vv(l)]; end ave_peroid=[ave_peroid sum(peroid)/(length(vv)-1)]; j1=j1+1; end HTT=[HTT abs(Nset-sum(ave_peroid)/N1)/Nset]; HTT1=[HTT1 abs(Nset-sum(ave_peroid1)/N1)/Nset]; %HTT2=[HTT2 abs(Nset-sum(ave_peroid2)/N1)/Nset]; end figure(12) semilogy(0:20,HTT,'b-',0:20,HTT1,'b*-') xlabel('信噪比/dB') ylabel('百分比误差') title('STFT周期估计误差') legend('优化算法','传统算法') toc; time=toc