%qpsk.m
%
%QPSK Simulation
%
%program by 07041061
%
function qpsk()
clear;
clc;
echo on
%*********************variable preparation******************************
ebn0=4;
sr=256000; %symbol rate
m1=2; %modulation lever
br=sr.*m1; %bit rate
nd =100; %number of symbols that simulates
%********************axes properity*************************************
T=1/sr; % symbol duration
delta_T=T/200; % sampling interval,compoversample function
alpha=0.5; % roll-off factor
fc=60/T; % carrier frequency
tt=10; % number of time period
t=-tt/2*T+delta_T:delta_T:tt/2*T; % time axis,N=2000
t1=-tt/2*(1+(nd-1)/tt)*T+delta_T:delta_T:tt/2*(1+(nd-1)/tt)*T;%t1的长度为109*200;每个时间点对应后面的扩充后的数组中的一个数;这个时间是自己设定的一个时间，因为在matlab里面没有时间的概念。只有离散的数点
N=length(t); % number of samples,N=2000
f1=-0.5*(1+(nd-1)/tt)/delta_T:1/(delta_T*N):0.5*(1+(nd-1)/tt)/delta_T-1/(delta_T*N);% actual frequency scale
f=-0.5/delta_T:1/(delta_T*(N-1)):0.5/delta_T;
% ******************calculation of raised cosine pulse*****************
for i=1:N,
if (abs(t(i))~=T/(2*alpha)),
g_T(i)=sinc(t(i)/T)*(cos(pi*alpha*t(i)/T)/(1-4*alpha^2*t(i)^2/T^2));
else
g_T(i)=0; % the value of g_T is 0 at t=T/(2*alpha)
end; % and at t=T(2*alpha)
echo off;
end;
echo on;
%********************fading initialnization******************************
tstp=1/sr/N;%
itau=[0];
dlvl=[0];
th1=[0.0];
now1=1;%只考虑了直射波多径所形成的瑞利衰落。
n0=[6];
itnd1=[1000];
fd=160;
flat=1;
%**************************data generration**************************
data=rand(1,nd*m1)>0.5;
disp(length(data));
%**********************qpsk source modulation**********************
paradata = data;
para =1;
m2=m1./2;
paradata2=paradata.*2-1;
count2=0;
for jj=1:nd,
isi = zeros(para,1);
isq = zeros(para,1);
for ii = 1:m2,
isi = isi + 2.^(m2 - ii).*paradata2((1:para),ii+count2);%串并变换I支路
isq = isq + 2.^(m2 - ii).*paradata2((1:para),m2+ii+count2);%串并变换Q支路
end;
iout((1:para),jj)=isi;
qout((1:para),jj)=isq;
count2=count2+m1;
end;
ich = iout;
qch = qout;
u_m_i=zeros(1,N*(1+(nd-1)/tt));% N/tt=T/delta_T=200;N=2000，扩大了为总长度为(1,109*200)的数组；
u_m_q=zeros(1,N*(1+(nd-1)/tt));%扩大数组的长度
for i=(1:nd)
u_m_i(1+(i-1)*N/tt:N+(i-1)*N/tt)=u_m_i(1+(i-1)*N/tt:N+(i-1)*N/tt)+ich(i)*g_T;% the modulated signal N/tt=T/delta_T=200;N=2000，这里200个数为一段
% g_T长度为2000(N) ，将每个符号的值与2000个成形函数的采样点相乘来做为结果的一段值，然后进行传输
u_m_q(1+(i-1)*N/tt:N+(i-1)*N/tt)=u_m_q(1+(i-1)*N/tt:N+(i-1)*N/tt)+qch(i)*g_T;% the modulated signal 进行脉冲成形为升余弦函数
end
%disp(length(u_m_i));
%*************************multiply carrier wave send*******************
u_m_i_cp1=u_m_i;u_m_q_cp1=u_m_q;%要画图的第一组数据
cos_carrier=cos(2*pi*fc*t1);%以上面自定的时间轴来产生载波离散点，长度为200*109
sin_carrier=sin(2*pi*fc*t1);%上下为两路正交载波
u_m_i=u_m_i.*cos_carrier;%两者的长度都为200*109;
u_m_q=u_m_q.*sin_carrier;
u_m_i_cp2=u_m_i;u_m_q_cp2=u_m_q;
%**********************Attenuation Calculation********************************************
spow=sum(ich.*ich+qch.*qch)/nd; %先计算信号强度
attn=0.5*spow*sr/br*10.^(-ebn0/10);%再根据信号强度来计算噪声幅度
attn=sqrt(attn);
%**********************Add White Gausion Noise(AWGN)********************************************
inoise=randn(1,length(u_m_i)).*attn;%noise
qnoise=randn(1,length(u_m_q)).*attn;%noise
u_m_iaw=u_m_i+inoise;%add awgn noise
u_m_qaw=u_m_q+inoise;%add awgn noise
% **********************fading******************************************************************
[ifade,qfade]=sefade(u_m_i,u_m_q,n0,itnd1,length(u_m_i),tstp,fd,flat);%itau,dlvl,th1,itnd1=1000,counter;
u_m_i=ifade+inoise;
u_m_q=qfade+qnoise;
%************************ receiver multiply carrier wave ******************************************
u_m_i_cp3=u_m_i;u_m_q_cp3=u_m_q;
u_m_i=u_m_i.*cos_carrier;%两者的长度都为200*109;
u_m_q=u_m_q.*sin_carrier;
u_m_i_cp4=u_m_i;u_m_q_cp4=u_m_q;
%*******************butterworth filter(cut-off frequency is chosen as same as carrier)****************
[b,a]=butter(21,fc/200/sr); %fc=60*sr;b,a只是返回的设计滤波器的两个参数,为一个行向量；
u_m_i=filter(b,a,u_m_i);%通过上面返回的两个行向量参数来设计滤波器再与待处理数据进行滤波运算；
u_m_q=filter(b,a,u_m_q);%同上；
u_m_i_cp5=u_m_i;u_m_q_cp5=u_m_q;
for i=1:nd
%I支路处理
for j=0:199
u_d(j+1)=u_m_i(N/2+(i-1)*N/tt+j);%分段处理结果数组
end;
index=find(abs(u_d)==max(abs(u_d)));%。找到这组数中绝对值最大数的索引号:index;对其进行保存以备下面的判决；
ich1(i)=u_d(index);
%Q支路处理
for j=0:199
u_d(j+1)=u_m_q(N/2+(i-1)*N/tt+j);%同上
end;
index=find(abs(u_d)==max(abs(u_d)));%同上
qch1(i)=u_d(index);
end
%************************ qpsk demodulation************************
idata = ich1;
qdata = qch1;
para = 1;
demodata=zeros(para,m1*nd);
demodata((1:para),(1:m1:m1*nd-1)) = idata((1:para),(1:nd))>=0;%I支路信号并串转换
demodata((1:para),(2:m1:m1*nd)) = qdata((1:para),(1:nd))>=0;%Q支路信号串并转换
% ************************ plotting commands follow ************************
figure(1);
subplot(2,1,1);
bar(data(1:10),'b');
title('Input baseband waveform');
DATA=abs(fft(data));
subplot(2,1,2);
plot(fftshift(DATA),'r');
title('power spectrum of Input baseband waveform');
figure(2);
subplot(2,1,1);
plot(t1,u_m_i_cp1);
title('插零后I支路滤波后信号,时域图');
subplot(2,1,2);
plot(f1,fftshift(abs(fft(u_m_i_cp1))));
title('插零后I支路滤波后信号,功率谱');
figure(3);
subplot(2,1,1);
plot(t1,u_m_i_cp2);
title('(QPSK modulated signal)加上载波后的I支路时域图，载频为基频的60倍');
subplot(2,1,2);
plot(f1,fftshift(abs(fft(u_m_i_cp2))));
title('(QPSK power spectrum)加上载波后的I支路功率谱图，载频为基频的60倍 ');
figure(41);
subplot(2,1,1);
plot(t1,u_m_iaw);%I支路信号出高斯噪声信道后的时域信号
title('AWGN channel output signal');
subplot(2,1,2);
plot(f1,fftshift(abs(fft(u_m_iaw))));%I支路信号出高斯噪声信道后的信号频谱
title('AWGN channel output power spectrum ');
figure(42);
subplot(2,1,1);
plot(t1,u_m_i_cp3);
title('(Rayleigh channel)接收到的I支路包含瑞利衰落加高斯噪声信号,时域图');
subplot(2,1,2);
plot(f1,fftshift(abs(fft(u_m_i_cp3))));
title('(Rayleigh channel)接收到的I支路包含瑞利衰落加高斯噪声信号，功率谱图');
figure(5);
subplot(2,1,1);
plot(t1,u_m_i_cp4);
title('接收到的I支路再次乘载波信号时域图');
subplot(2,1,2);
plot(f1,fftshift(abs(fft(u_m_i_cp4))));
title('接收到的I支路再次乘载波信号功率谱');
figure(6);
subplot(2,1,1);
plot(t1,u_m_i_cp5);
title('接收到的I支路滤波后信号,时域图');
subplot(2,1,2);
plot(f1,fftshift(abs(fft(u_m_i_cp5))));
title('接收到的I支路滤波后信号,功率谱');
figure(7);
subplot(2,1,1);
bar(demodata(1:10),'b');
title('Demodulated signal');
DEMODATA=abs(fft(demodata));
subplot(2,1,2);
plot(fftshift(DEMODATA),'r');
title('Power spectrum of Demodulated signal');
scatterplot([ich' qch']);
title('输入信道前的星座图 ');
scatterplot([u_m_i_cp1' u_m_q_cp1']);
title('插零以及脉冲成形后的星座图 ');
scatterplot([u_m_i_cp2' u_m_q_cp2']);
title('进入瑞利信道前的星座图 ');
scatterplot([u_m_i_cp3' u_m_q_cp3']);
title('出瑞利信道后的星座图 ');
scatterplot([ich1' qch1']);
title('解调后的星座图 ');
%************************ end ************************
评论0