psk_cyclo_corr.m

psk 的编码调试,系统性能测试,利用MATLAB编程实现,仅供参考

%%%  设置调制参数  bpsk qpsk msk
fc=200;
fd=40;
fs=1000;

%%%  生产带通信号

% s=randint(1,400,[0,1]);
% y=dmod(s,fc,fd,fs,'psk',2);
% h1=fir1(50,[(fc-fd)/fs*2,(fc+fd)/fs*2]);
% y=conv(y,h1);

% s=randint(1,400,[0,3]);
% y=dmod(s,fc,fd,fs,'psk',4);
% h1=fir1(50,[(fc-fd)/fs*2,(fc+fd)/fs*2]);
% y=conv(y,h1);

% s=randint(1,400,[0,1]);
% y=dmod(s,fc,fd,fs,'msk');

%%% 设置调制参数 am fm
fc=200;
fs=600;
%%%  生产带通信号
n=randn(1,2200);
Xm(1)=n(1);
for i=2:2200
    Xm(i)=0.979274*Xm(i-1)+n(i);
end
y=amod(Xm,fc,fs,'amdsb-sc');
% y=amod(Xm,fc,fs,'fm');

%%%  生产解析信号
h=hilbert(y);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  产生用于dsp调试的数据
N=length(h);
fid=fopen('input_data.dat','w');
fprintf(fid,'1651 4 0 0 0\r\n');
for i=1:N
    fprintf(fid,'%f\r\n',real(h(i)));
    fprintf(fid,'%f\r\n',imag(h(i)));
end
fclose(fid);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  计算循环自相关  2048点fft
spec=[];
J=zeros(1,2048);
M=10;
x=h(129:129+2047);
for tao=0:M-1
  xn=h(129-tao:129-tao+2047);
  z=x.*conj(xn);           % 共轭相乘
  spec1=abs(fft(z)).^2;    % 计算频谱
  J=J+spec1;               % 频谱累加
  spec=[spec;spec1];       % 存储频谱，用于调试时观测
end

%%% 画循环自相关
subplot(2,1,1);
mesh(spec(:,2:1024));      % 为了观测方便，去掉第一列（对应的循环频率为0）
subplot(2,1,2);
plot(J(2:1024));           % 为了观测方便，去掉第一列（对应的循环频率为0）

%%% 识别
thro1=1;                   % 判决门限
JJ=J(1:1024);JJ(1)=0;      % 找到2--1024之间最高的谱线
[val alpha]=max(JJ) ;      
                           

% %%% 第一种判决规则
% N=5;                       % 取最高谱线左右两边各N根谱线之和做归一化
% t=[alpha,alpha+1:alpha+N];       % 找到用于归一化的下标
% for i=1:N
%     if alpha-i==0
%         temp=2048;
%     else
%         temp=mod(alpha-i,2048);
%     end
%     t=[temp t];
% end
% num=sum(J(t(5:7)));
% den=sum(J(t))-num;
% Rx=num/den                 % 归一化
% esti_fd=(alpha-1)*fs/2048;
% 
% if Rx < thro1
%     str='it is am,fm or msk' 
% else
%     str='it is bpsk or qpsk'
%     esti_fd=(alpha-1)*fs/2048
% end


%%% 第二种判决规则
if alpha<10+1+1
%     str='it is am,fm or msk'
    env=abs(h(1:2048)).^2;
    m_env=mean(env);
    env_cn=(env/m_env)-1;
    env_fft=abs(fft(env_cn)).^2;    % 此处不应该平方，应该除以2048，为了降低dsp计算复杂度，省略此除法运算
    env_fft(1)=0;
    max_env=max(env_fft);
    if max_env<1000
        str='it is fm or msk'
    else
        str='it is am'
    end
else
    t=[alpha-10:alpha+10];  % 生成用于归一化的下标
    num=sum(J(t(10:12)));
    den=sum(J(t))-num;
    Rx=num/den;                 % 归一化
    if Rx < thro1
%         str='it is am,fm or msk' 
        env=abs(h(1:2048)).^2;
        m_env=mean(env);
        env_cn=(env/m_env)-1;
        env_fft=abs(fft(env_cn)).^2;    % 此处不应该平方，应该除以2048，为了降低dsp计算复杂度，省略此除法运算
        env_fft(1)=0;
        max_env=max(env_fft);
        if max_env<1000
            str='it is fm or msk'
        else
            str='it is am'
        end
    else
        esti_fd=(alpha-1)*fs/2048
        h2fft=abs(fft(h(1:2048).^2));         % 平方，FFT
        m_h2fft=mean(h2fft);
        h2fft_cn=(h2fft/m_h2fft)-1;
        [max_h2fft index]=max(h2fft);
        num=sum(h2fft(index-2:index+2));
        den=sum(h2fft(index-10:index+10))-num;
        rfc=num/den;
        if rfc>1
            esti_fc=(index-1)/2048*fs/2
            str='it is bpsk'
        else
            h4fft=abs(fft(h(1:2048).^4));         % 平方，FFT
            [max_h4fft index]=max(h4fft);
            esti_fc=(index-1)/2048*fs/4
            str='it is qpsk'
        end
    end
end



h2fft=abs(fft(h(1:2048)));         % 平方，FFT
plot(h2fft)
% m_h2fft=mean(h2fft);
% h2fft_cn=(h2fft/m_h2fft)-1;
% [max_h2fft index]=max(h2fft);
% num=sum(h2fft(index-2:index+2));
% den=sum(h2fft(index-10:index+10))-num;
% rfc=num/den
% 
% % var(h2fft_cn)
% plot(h2fft_cn)