mfsk5.m

对数字调制信号计算循环谱相关

% =============================================
% ----------------BP_MFSK_simulation----------
% --------------频移键控信号仿真程序
% --------------程序内部可设定码速率等参数
% =============================================
clear all;

Fc=1.0e+8;%载频  
Fsc=4*Fc;
Fd=2.5e+7;%码元速率
Fs=4.0e+8;       %--------系统采样频率  !!!!!!!!
tone=5.0e+7 ;
    %---------MFSK-------
    m=2;
numSymb=32  ;%进行仿真的信息代码个数
x=randsrc(numSymb,1,[0:m-1]);%---------产生25个二进制随机码
x1=dmod(x,Fc,Fd,Fs,'fsk',m,tone);
x2=dmod(x,Fc,Fd,Fs,'psk',m);
x3=dmod(x,Fc,Fd,Fs,'ask',m);
x4=dmod(x,Fc,Fd,Fs,'msk');

numModPlot=numSymb*Fs/Fd;%－－－－－－注意调制的点数－－－－－－－－－
t=[0:numModPlot-1]'./Fs;
figure(1)
subplot(411);
plot(t,x1(1:length(t)),'b-');
axis([min(t) max(t) -1.5 1.5]);
title('fsk调制后的信号');
xlabel('Time');
ylabel('Amplitude');

subplot(412);
plot(t,x2(1:length(t)),'b-');
axis([min(t) max(t) -1.5 1.5]);
title('psk调制后的信号');
xlabel('Time');
ylabel('Amplitude');

subplot(413);
plot(t,x3(1:length(t)),'b-');
axis([min(t) max(t) -1.5 1.5]);
title('ask调制后的信号');
xlabel('Time');
ylabel('Amplitude');

subplot(414);
plot(t,x4(1:length(t)),'b-');
axis([min(t) max(t) -1.5 1.5]);
title('msk调制后的信号');
xlabel('Time');
ylabel('Amplitude');

figure(2)

n1=length(x1);
y1=fft(x1)/n1;
f=(-Fs/2:Fs/n1:Fs/2-Fs/n1);
subplot(4,1,1);
plot(f,fftshift(abs(y1)));
title('FSK信号FFT频谱');
xlabel('Frequency');
ylabel('Amplitude(归一化后)');

n2=length(x2);
y2=fft(x2)/n2;
f=(-Fs/2:Fs/n2:Fs/2-Fs/n2);
subplot(4,1,2);
plot(f,fftshift(abs(y2)));
title('PSK信号FFT频谱');
xlabel('Frequency');
ylabel('Amplitude(归一化后)');

n3=length(x3);
y3=fft(x3)/n3;
f=(-Fs/2:Fs/n3:Fs/2-Fs/n3);
subplot(4,1,3);
plot(f,fftshift(abs(y3)));
title('ASK信号FFT频谱');
xlabel('Frequency');
ylabel('Amplitude(归一化后)');


n4=length(x4);
y4=fft(x4)/n4;
f=(-Fs/2:Fs/n4:Fs/2-Fs/n4);
subplot(4,1,4);
plot(f,fftshift(abs(y4)));
title('MSK信号FFT频谱');
xlabel('Frequency');
ylabel('Amplitude(归一化后)');

%----------------------加高斯白噪声---------------------
figure(3)
SNRperbit=5;
adjSNR=SNRperbit-10*log10(Fs/Fd)+10*log10(log2(m));

Y1=awgn(x1,adjSNR, 'measured', [], 'dB');
subplot(411);
plot(t,Y1(1:length(t)),'r');
title('FSK信号加高斯白噪声信号');
ylabel('Amplitude');
xlabel('Time');

Y2=awgn(x2,adjSNR, 'measured', [], 'dB');
subplot(412);
plot(t,Y2(1:length(t)),'r');
title('PSK信号加高斯白噪声信号');
ylabel('Amplitude');
xlabel('Time');

Y3=awgn(x3,adjSNR, 'measured', [], 'dB');
subplot(413);
plot(t,Y3(1:length(t)),'r');
title('ASK信号加高斯白噪声信号');
ylabel('Amplitude');
xlabel('Time');

Y4=awgn(x4,adjSNR, 'measured', [], 'dB');
subplot(414);
plot(t,Y4(1:length(t)),'r');
title('MSK信号加高斯白噪声信号');
ylabel('Amplitude');
xlabel('Time');

figure(4)
N1=fft(x1)/n1;
subplot(4,1,1);
plot(f,fftshift(abs(N1)));
title('FSK信号加高斯白噪声信号FFT频谱');
xlabel('Frequency');
ylabel('Amplitude(归一化后)');

N2=fft(x2)/n2;
subplot(4,1,2);
plot(f,fftshift(abs(N2)));
title('PSK信号加高斯白噪声信号FFT频谱');
xlabel('Frequency');
ylabel('Amplitude(归一化后)');

N3=fft(x3)/n3;
subplot(4,1,3);
plot(f,fftshift(abs(N3)));
title('ASK信号加高斯白噪声信号FFT频谱');
xlabel('Frequency');
ylabel('Amplitude(归一化后)');


N4=fft(x4)/n4;
subplot(4,1,4);
plot(f,fftshift(abs(N4)));
title('MSK信号加高斯白噪声信号FFT频谱');
xlabel('Frequency');
ylabel('Amplitude(归一化后)');


%-----------------------------------------------------
%周期图时域平滑法求循环谱相关


%figure(2)
%N=numModPlot;M=64;s1=y';
%S=cyclic_cross_periodogram1(s1,s1,N,M,'hamming','hamming');
%-----------------------------------------------
%周期图频域平滑法求循环谱相关
N=2048;
M=4;
K=round(1.5*N*Fc/Fsc);          %频率轴
alfa=round(1.5*N*Fc/Fsc);  %循环频率轴 

G1=cyclic_function(Y1,N,M,K,alfa,Fs,Fsc);
G2=cyclic_function(x2,N,N,K,alfa,Fs,Fsc);
G3=cyclic_function(x3,N,M,K,alfa,Fs,Fsc);
G4=cyclic_function(x4,N,N,K,alfa,Fs,Fsc);



%X=fft(s2,N);
%for k=1:5*N
 %  if k>=1&k<=N
  %    XX(k)=X(k);
  % else if mod(k,N)~=0
   %        XX(k)=X(k-floor(k/N)*N);
    %   else XX(k)=X(N);
    %   end
   %end
%end

 
%K=round(1.5*N*Fc/Fsc);          %频率轴
%alfa=round(1.5*N*Fc/Fsc);       %循环频率轴  
% K=round(2*N*fc/fsc);
% alfa=round(4*N*fc/fsc);
%i=1;
%for k=-K:16:K
    %j=1;
    %for m=-alfa:16:alfa
     %   sum=0;
     %   for v=-M:M
      %      a=k+v+m;b=k+v-m;
       %     if (k+v+m>=0)&&(k+v-m>=0)
        %        sum=sum+XX(k+v+m+1)*conj(XX(k+v-m+1));
         %   else while(a<0||b<0)
          %      a=a+N;
           %     b=b+N;
            %    end
             %   sum=sum+XX(a+1)*conj(XX(b+1));
            %end
        %end
        
        %g_en(i,j)=sum/Fs/2/M/N;
      
        %j=j+1;
    %end
    %i=i+1;
%end

%p=(-K:16:K)*Fsc/N;
%q=(-alfa:16:alfa)*2*Fsc/N;

%p=[-0.5:1/N:0.5-1/N]*2*pi;
%q=(-0.5:1/N:0.5-1/N)*2;


%figure(2)
%G=g_en.*g_en;
%z2=zeros(length(q));
%mesh(p,z2,abs(G));
%Sm=max(max(abs(G)));


%figure(2)
%mesh(p,q,abs(g_en));

%figure(3)
%subplot(2,2,2);
  %plot(q*N/2,abs(g_en));
%z2=zeros(length(p));
%mesh(z2,q,abs(g_en));








%--------------------------------------------------------------------------
%---在0频率轴上对循环频率轴进行一维搜索估计载频-------------------------------
%i=1;
%while((-K+(i-1)*32)~=0)
 %   i=i+1;
%end

%for j=1:length(-alfa:32:alfa)
 %   f0(j)=g_en(i,j);
%end

 %[mm,ff]=max(f0);
 %(-alfa+(ff-1)*32)*2*Fsc/N
 
 %-------------------------------------------------------------------------