qpsk_awgn_draw.m

matlab的调制解调

% 这个仿真程序要重点掌握：
% 其他的程序没有考虑载波的问题；

% 总结一下 QPSK 的仿真过程：产生基带信号，一个符号有很多的采样
% 点数这个点数要设好，编码，串并变换，产生载波信号，
% 载波信号的采样间隔和基带符号的采样间隔要一样，上述两者相乘
% 就完成了调制，加噪声，接受信号和发送端载波相乘就完成解调，
% 判决，并串变换，统计并计算误码率； 


% 这里要尤为注意一点就是：基带信号的采样间隔和载波信号的采样间隔
% 要设为一样的，这样便于观察频谱的搬移；

% QPSK with AWGN channel 
clear all; close all;
% initialization
T=1;                          % symbol duration
fc=10/T;                      % carrier frequence

% modulation of level ml＝1 BPSK，ml＝2 QPSK
ml=2;
% number of transmiting bit，发送的比特数；
nb=100;                       
% sampling interval，抽样间隔，一个符号抽样200个点;
% 对一个符号抽样多个点，这是以前常常忽视的。
% 这里基带信号的采样频率为200Hz，可以看到后面的载波信号的
% 采样频率也是200Hz；
delta_T=T/200;                
% sampling frequency，抽样频率； 
fs=1/delta_T;                 
SNR=0;                        % signal to noise ratio
t=0:delta_T:nb*T-delta_T;
N=length(t);                  % number of samples  
% generate source data
data=randn(1,nb)>0.5;  
% translate source data into NRZ code
% 注意是对符号编码；
datanrz=data.*2-1;         
% 虽然发送了100个符号，但每个符号的采样点数为200，共有20000个点；
data1=zeros(1,nb/delta_T); 
% data1是基带发送信号；
for q=1:nb
    data1((q-1)/delta_T+1:q/delta_T)=datanrz(q);
end;                           % transmitting signal 

% 是基带发送信号的FFT变换；
data2=abs(fft(data1));

% S/P conversion
% 串并变换；
idata=datanrz(1:ml:(nb-1));
qdata=datanrz(2:ml:nb);

% QPSK modulation
ich=zeros(1,nb/delta_T/2);
for i=1:nb/2
    ich((i-1)/delta_T+1:i/delta_T)=idata(i);
end;
% 符号的持续时间是T，要使在符号的持续时间内调制后的信号功率为1
% 必须对载波加 sqrt(2/T) 的幅度；
% 以下是同相信道调制；
for ii=1:N/2
    a(ii)=sqrt(2/T)*cos(2*pi*fc*t(ii));
end;
idata1=ich.*a;                % I-channel multiply carrier
% 以下是正交信道调制；
qch=zeros(1,nb/2/delta_T);
for j1=1:nb/2
    qch((j1-1)/delta_T+1:j1/delta_T)=qdata(j1);
end;
for jj=1:N/2
    b(jj)=sqrt(2/T)*sin(2*pi*fc*t(jj));
end;
qdata1=qch.*b;               % Q-channel multiply carrier

% QPSK modulated signal
s=idata1+i*qdata1;             
ss=abs(fft(s));

% 在这个信道中只考虑了加性高斯噪声；
% through AWGN 
s1=awgn(s,SNR);              % received signal
s11=abs(fft(s1));

% QPSK demodulation
% QPSK的解调部分，就是将接受信号和发送端的载波信号相乘；
idata2=s1.*a;
qdata2=s1.*b;                % multiply carrier

idata3=zeros(1,nb/2);
qdata3=zeros(1,nb/2);
% 以下是解调后的判决部分；
for n=1:nb/2
    if sum(idata2((n-1)/delta_T+1:n/delta_T))>=0
        idata3(n)=1;
    else idata3(n)=0;
    end;
    if sum(qdata2((n-1)/delta_T+1:n/delta_T))>=0
        qdata3(n)=1;
    else qdata3(n)=0;
    end
end                         % integral and decision

% P/S conversion
% 再进行并串变换；
demodata=zeros(1,nb);
demodata(1:ml:(nb-1))=idata3;
demodata(2:ml:nb)=qdata3;

% calculate number of error
% 统计误码个数；
num_BER=sum(abs(demodata-data))

% plot commends
h = spectrum.welch;

figure(1);
plot(data1);
title('baseband signal');

figure(2);
plot(fs*(-10000:10000-1)/20000,fftshift(data2));
title('baseband signal spectrum');

figure(3);
psd(h,data1,'fs',fs);
title('baseband signal power spectrum');

figure(4);
plot(s);
title('modulated signal');
axis([0 1000 -3 3]);

figure(5);
plot(fs*(-5000:5000-1)/10000,fftshift(ss));
title('modulated signal spectrum');

figure(6);
psd(h,s,'fs',fs);
title('modulated signal power spectrum');

figure(7);
plot(s1);
title('AWGN output signal');
axis([0 500 -6 6]);

figure(8);
plot(fs*(-5000:5000-1)/10000,fftshift(s11));
title('output signal spectrum');

figure(9);
psd(h,s1,'fs',fs);
title('output siganl power spectrum');