qpsk_rayleigh_draw.m

matlab的调制解调

% 本程序相对于程序QPSK_AWGN_draw加入了瑞利衰落，其他的都是一样的
% 缺点：没有考虑多径，频移；
% QPSK with Rayleigh channel
clear all;; close all;
% initialization
% symbol duration
T=1;
% signal to noise ratio
SNR=0;
% carrier frequence
fc=6/T;
% modulation of level ml＝1 BPSK，ml＝2 QPSK
ml=2;
% number of transmiting bits
nb=100;
% sampling interval
delta_T=T/100;
% sampling frequency 
% 注意采样频率是 100 Hz；
fs=1/delta_T;
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;
data1=zeros(1,nb/delta_T);
for q=1:nb
    data1((q-1)/delta_T+1:q/delta_T)=datanrz(q);
end;                         % transmiting signal 
data2=abs(fft(data1));
% S/P conversion    
idata=datanrz(1:ml:(nb-1));
qdata=datanrz(2:ml:nb);

% QPSK modulation
ich=zeros(1,nb/2/delta_T);
for i=1:nb/2
    ich((i-1)/delta_T+1:i/delta_T)=idata(i);
end;
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 j=1:nb/2
    qch((j-1)/delta_T+1:j/delta_T)=qdata(j);
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

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

% through rayleigh channel
ray_ich=raylrnd(0.8,1,nb/2/delta_T);
ray_qch=raylrnd(0.8,1,nb/2/delta_T);
idata2=idata1.*ray_ich;
qdata2=qdata1.*ray_qch;
s1=idata2+qdata2;
s2=awgn(s1,SNR,'measured');

% QPSK demodulation
% 解调：
idata3=s2.*a;
qdata3=s2.*b;                 % multiply carrier

idata4=zeros(1,nb/2);
qdata4=zeros(1,nb/2);
for n=1:nb/2
    ichsum(n)=sum(idata3((n-1)/delta_T+1:n/delta_T))*delta_T;
    if ichsum(n)>=0
        idata4(n)=1;
    else idata4(n)=0;
    end
    qchsum(n)=sum(qdata3((n-1)/delta_T+1:n/delta_T))*delta_T;
    if qchsum(n)>=0
        qdata4(n)=1;
    else qdata4(n)=0;
    end
end                          % integral and decision

% P/S conversion
demodata=zeros(1,nb);
demodata(1:ml:(nb-1))=idata4;
demodata(2:ml:nb)=qdata4;

% 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*(-5000:5000-1)/10000,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 500 -3 3]);

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

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

figure(7)
for i=1:nb/2
plot(idata(i),qdata(i),'ro');
title('constellation');        hold on;
axis([-2 2 -2 2]);
plot(ichsum(i),qchsum(i),'*');
hold on;
legend('before channel','after channel');
end;