PAM信号的产生及功率谱计算（维纳-辛钦定理）

clc;
clear all
%%产生PAM
clear all;close all;clc;
n_symbol=5;     %%产生矩形脉冲个数
N = 100;   %每个符号的采样点数
T=1;                                %%设置时间为1，即矩形PAM的总功率应该为1
t=0:T/N:T-T/N;
num_cos = 4;    %确定每个脉冲用几个正弦表达
fc = num_cos/2; %因为频率要用全波                    
c=sqrt(2/T)*sin(pi*num_cos*t);      %载波信号
c = abs(c);                      %半波正弦
M=2;       %电平数量
graycode=[0 1 3 2];    %设置编码规则
noise = 0;   %是否设置噪声
msg=randi(M,1,n_symbol);
msg1=graycode(msg);                 %%按照之前定义的graycode规则重新定义各数
msgmod=pammod(msg1,M).';
x=msgmod*c;
pamsig=reshape(x.',1,length(msgmod)*length(c));
if noise == 1
    s_pow=norm(pamsig).^2/(n_symbol);     %%这里是符号的平均功率，而不是每一个采样点
    sigma=sqrt(s_pow/(2*snr));         %%根据线性值的信噪比和信号功率求出噪声功率
    pamsig=pamsig+sigma*randn(1,length(pamsig));              %%加性噪声后的信号
end

%画图拓宽时域
for i=1:n_symbol                    
    orisig((N*i-N+1):N*i)=msg1(i);
end
figure(1);
subplot(311);
plot(orisig);title('矩形脉冲调制PAM信号');
subplot(312);
plot(pamsig);title('正弦半波调制PAM信号');

figure(2)
not = 1024;
subplot(2,1,1)
[c,d] = xcorr(orisig); %%矩形脉冲方式的自相关
amp = fftshift(abs(fft(c,not)));
amp = amp/not;         %这里不是除以N/2，原因是这里我得到的是有负频率的谱
amp = amp/length(orisig);
plot((-not/2:not/2-1), smooth(20*log10(amp)));
title('矩形脉冲调制PAM信号功率谱')
subplot(2,1,2)
[a,b] = xcorr(pamsig); %%正弦方式自相关函数
ampp = fftshift(abs(fft(a,not)));
ampp = ampp/not;
ampp = ampp/length(pamsig);
plot(fc+(-not/2:not/2-1), smooth(20*log10(ampp))); %%平滑处理
title('正弦半波调制PAM信号')

po = sum(orisig.*orisig/length(orisig));
pog = sum(amp);
pp = sum(pamsig.*pamsig/length(pamsig));
ppg = sum(ampp);

%%功率谱正确性验证
%注意pamsig，因为fft只做了1024个点，n_symbol最多取5，大于5的话，要变大not
if round(po/pog) == 1
    disp('矩形脉冲PAM功率谱计算正确，功率为：');
    disp(pog);
end

if round(pp/ppg) == 1
    disp('半波正弦脉冲PAM功率谱计算正确，功率为：');
    disp(ppg);
end