est_symbol.m

基于小波变换的PSK信号符号率估计的matlab程序

clc;clear;

M          = 1200;               %%采样点数
loop       = 1000;                 %%循环次数
fs         = 10e6;               %%采样率
fc         = 2.5e6;              %%载频 
fcode      = 1e6;                %%符号率
fm         = fcode;              
N          = 1500;
SNR_array  = -5:15;              %%信号的信噪比


for SNR_index=1:length(SNR_array)
        SNR=SNR_array(SNR_index)-7;
        fprintf('processing... %d/%d times\n',SNR_index,length(SNR_array));
        for loop_index=1:loop    
            
            sig =sigmpsk(fc,fcode,fs,2,N);  % 产生一个信号, 并hilbert变换成复数形式
            sig=awgn(sig,SNR,'measured');
            sig=hilbert(sig);
            
            frqsig=fc/fs;
            morfrq=centfrq('morl');
            chdyz=morfrq/frqsig;                              %尺度因子值    
            
            b1=abs(cwt(sig,chdyz,'morl'));   %进行小波变换
            b1=b1-mean(b1); 
            
            c1=autocorr(b1,M+100); 
            c1=diff(diff(c1));              %对小波系数进行自相关，且进行二次差分  
                
            c1=c1(1:M);                     %对自相关值估计符号率
            H=hann(length(c1));
            s=hilbert(c1);
            p=fft(s.*H',M);
            B=abs(p);         %fft输出幅值%
                
            B=B(1:(M-10)/2);    %频谱中幅度最高的三根谱线及对应fft样点序号%
            F=find(B==max(B(2:end)));   %最高谱线%
            if not(rem(F-1,2))
               if B((F-1)/2+1)>0.6*B(F)
                  F=(F-1)/2+1;
               else
                  F=F;
               end
            else 
                F=F;
            end
 
            G=F-1;I=F+1;
            df1=-B(G)/(B(G)+B(F))*(F-G)*fs/M;    %基于两点插值的频率修正值计算%
            df2=B(I)/(B(I)+B(F))*(I-F)*fs/M;
            k1=(B(F)-B(G))/((F-G));              %相邻谱线幅值变化斜率%
            k2=(B(I)-B(F))/((I-F));
            W1=abs(k2)/(abs(k1)+abs(k2));        %三点插值的权系数计算%
            W2=abs(k1)/(abs(k1)+abs(k2));
            frq=(F-1)*fs/M+W1*df1+W2*df2;        %三点插值频率估计%
            
            err_codef(loop_index,SNR_index)=(frq-fcode)/fcode   %计算符号率估计相对误差
        end
end

%估计符号率估计均方根误差
for SNR_index=1:length(SNR_array)
err_rmse(SNR_index)=sqrt(mean(err_codef(:,SNR_index).^2));
end

%对统计结果绘图
figure;semilogy(SNR_array,err_rmse);title('符号率估计误差随SNR关系图(morl小波,1000次)')
  xlabel('信噪比/dB');ylabel('符号率归一化均方根误差');grid on;