基于Matlab实现跳频信号主要参数估计，包括跳变周期，跳时刻，瞬时频率估计，以及误差分析.zip

%输入信息序列，4FSK调制1111111111111111 x=input('Please enter :');%输入任意双比特数组，如：[11 01 10 01 00 01 10 ...]，但数组长度不大于10 N=length(x); Ts=0.00001; fs=1/Ts; Tm=0.002; EndTime=0.004-Ts; %2s Nset=(EndTime+Ts)/Ts; y=zeros(Nset,1); i=0; xd=zeros(1,N); 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]; b=input('选择跳频序列：');%输入1~9任意数字，选择不同跳频序列（最好选择4、5、6，这样保证宽跳频，每一跳频率差更大） Su=mod(b.*G,11); for j=2:11 fprintf('%3d',Su(j)) end % 跳频序列合成跳频载波1111111111111111111 Tc=0.0003; C=zeros(Nset,1); for k=2:(N+1) [Carrier1,time2]=gensig('sin', 1/(Su(k)/Tc+xd(k-1)/Tm),EndTime,Ts); C(1+(k-2)*Nset:Nset*(k-1),1)=Carrier1; end %混频产生跳频信号 FHsignal=C; subplot(2,1,2) plot(FHsignal) xlabel('t') ylabel('跳频调制信号') title('跳频信号') grid %1111111111111111111111111111111111111111111111 z=hilbert(C); e=2^nextpow2(length(z)/4); h1= window('hamming',401); [tfr ,t,f]= TFRSTFT(z,1:length(z),e,h1,0); %%修改：1、我的短时傅里叶变换函数是tfrstft，你在你那里运行时改回TFRSTFT；2、短时傅里叶变换出来的频率有500个点 [fmax,ff]=max(abs(tfr)); fFH1=abs(fft(fmax)); [maxr1,f1]=max(fFH1(2:length(fFH1))); Nlag=length(z)/f1; Fo=abs(fft(z)); width=50; high=max(Fo)*0.8; for M=1:length(Fo) if (Fo(M)<high) Fo(M)=0; end end [peak,loc] = findpeaks(Fo,'minpeakdistance',width); Lloc=length(loc); loc1=loc(2:Lloc)-loc(1:Lloc-1); N1=min(loc1); Ndop=ceil(4000*1.5/N1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % apauto=autocorr(z,length(z)-1); % % % auto=abs(apauto); % % % figure(2) % % % plot(auto); % % % d=0; % % % sumd=zeros(1,100); % % % while(d<100) % % % sumd(d+1)=auto(1+10*d)+auto(2+10*d)+auto(3+10*d)+auto(4+10*d)+auto(5+d*10); % % % sumd(d+1)=auto(6+10*d)+auto(7+10*d)+auto(8+10*d)+auto(9+10*d)+auto(10+10*d); % % % d=d+1; % % % end % % % [ad,ld]=min(sumd); % % % Nlag=ld*10+1; % % % Nd(cc+7)=Nlag; % % % apf=fft(z); % % % figure(3) % % % plot(abs(apf)) % % % % % % % % end % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 11111111111111111111111111111111111111111111111111111111111111111111 % 计算SWWVP分布1111111111111111111111111 if(mod(Nlag,2)==0) Nlag=Nlag+1; end if(mod(Ndop,2)==0) Ndop=Ndop-1; end h2= window('hanning',Ndop); h3= window('hamming',Nlag); e=2^nextpow2(length(z)/4); [tfr2,t2,f2]=TFRSPWV(z,1:length(z),e/4,h2,h3,0); figure(4) contour(t2,f2,abs(tfr2)); set(gca,'ydir','normal') xlabel('shijian t'); ylabel('归一化频率'); ylim([0 0.5]) title('SWWVD时频分析图') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%参数估计：瞬时频率法%%%%%%%%%%%%%%%%%%%% [Af,locf]=max(tfr2); figure(5) plot(t2,locf/512) xlabel('时间T') ylabel('归一化频率F') ylim([0 0.5]) title('SWWVD时频脊线') locf1=locf(201:length(locf)); locf2=locf1(1:length(locf1)-1)-locf1(2:length(locf1)); locf22=abs(locf2); paixu=sort(locf22,'descend'); Ti=find(locf22>paixu(10))+200;%%%%%%%%%%%%%%%%%%跳时刻估计 LTi=length(Ti); Dt=(Ti(2:LTi)-Ti(1:LTi-1))*ones(LTi-1,1)/(LTi-1);%%%%%%%%%%%%%%%跳周期估计 locfbu=[locf,0,256]; Fi=hist(locfbu,e/4); Fmax=max(Fi)*0.8; Fgui=find(Fi>Fmax)*2/e;%%%%%%%%%%%%%%%%%%%%跳频率估计：归一化频率集 Ffi=Fgui/Ts; %%%%%%%%%%%%%%%%%%%%%%%%转换后的真实频率集 NLf=length(Ffi); ai=1; Ti1=[1,Ti,length(locf)]; sum1=zeros(1,NLf); IF=zeros(1,length(locf)); ccc=zeros(1,20); while(ai<NLf+1) dt=Ti1(ai+1)-Ti1(ai); sum1(ai)=locf(Ti1(ai):Ti1(ai+1)-1)*ones(dt,1)/dt*2/e;%%%%%%%%%归一化平均频率 [mi,locmi]=min(abs(Fgui-sum1(ai))); ccc(ai)=locmi; IF(Ti1(ai):Ti1(ai+1)-1)=Fgui(locmi);%%%%%%%%%%%瞬时频率估计 ai=ai+1; end figure(6)%%%%%%%%%%%% plot(IF) xlabel('时间T') ylabel('归一化频率F') ylim([0 0.5]) title('SWWVD瞬时频率估计') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% result=zeros(9,9); err0=201:400:3401; err1=zeros(9,9); HT=zeros(1,9); aba=1; for i=-8:0 j=0; vv1=zeros(1,9); bei=0; while(j<10) y3 = awgn(C,i,'measured'); z3=hilbert(y3); %%%%%%%%%%% [tfr3 ,~,~]= TFRSTFT(z3,1:length(z3),e,h1,0); %%修改：1、我的短时傅里叶变换函数是tfrstft，你在你那里运行时改回TFRSTFT；2、短时傅里叶变换出来的频率有500个点 [fmax,~]=max(abs(tfr3)); fFH1=abs(fft(fmax)); [~,f1]=max(fFH1(2:length(fFH1))); Nlag=length(z3)/f1; if(mod(Nlag,2)==0) Nlag=Nlag+1; end h4= window('hanning',101); h5= window('hamming',Nlag); [tfr4,t,f]=TFRSPWV(z3,1:length(z),e/4,h4,h5,0); %%%%%%%%%%%%%%%%%%%%% [fmax,ff]=max(abs(tfr4)); fFH1=abs(fft(fmax)); [maxr1,f1]=max(fFH1(2:length(fFH1))); EN1=length(z)/f1; ht=ff(2:length(ff)); ht1=ff(1:length(ff)-1)-ht; htt2=abs(ht1(201:length(ht1)-200)); a=max(htt2)*0.3; ah=htt2; for M=1:length(htt2) if (htt2(M)<a) htt2(M)=0; end end [~,vv] = findpeaks(htt2,'minpeakdistance',200); if(length(vv)==9) vv1=abs(vv-err0)+vv1; HT(i+9)=abs((vv(2:9)-vv(1:8))*ones(8,1)/8-400)+HT(i+9);%%%%%%tiaozhouqiwucha else bei=bei+1; end j=j+1; end vv2=vv1/(10-bei);%%%%%%%%%tiaoshke HT(i+9)=HT(i+9)/(10-bei); result(:,i+9)=vv2; aba=aba+1; end result1=result/4; figure(8) HT1=HT/4; plot(-8:0,HT1) xlabel('信噪比db') ylabel('百分比误差') title('SWWVD周期估计误差') figure(9) plot(abs(ht1)) xlabel('时间t'); ylabel('幅度'); title('SWWVD跳时刻') figure(10) mesh(-8:0,1:9,result1) xlabel('信噪比db'); ylabel('相对跳时刻'); title('SWWVD跳时刻估计误差百分比') aver=mean(result1); figure(11) plot(-8:0,aver,'-r') xlabel('信噪比db') title('SWWVD跳时刻平均误差') figure(12) a9=htt2(vv); stem(vv,a9) xlabel('时间T') title('SWWVD跳时刻估计值')