%% QPSK
clc;clear all;close all;
%假定接收端已经实现载波同步,位同步
%(盲信号解调重点要解决的问题:载波同步(costas环(未见到相关代码)),位同步(Gardner算法(未见相关代码)),帧同步)
%carrier frequency for modulation and demodulation
fc=5e6; %码数率为5MHZ
%QPSK transmitter
data=100; %原码个数
%随机数产生单极性码,randn返回1*100的数组,标准正态分布的标量
rand_data=randn(1,data);
for i=1:data
if rand_data(i)>=0.5
rand_data(i)=1;
else
rand_data(i)=0;
end
end
num_fig = 1;
figure(num_fig);%创建窗口
num_fig = num_fig+1;
subplot(211)
%画出随机序列的图像
stem(rand_data);%将数据序列 rand_data 绘制为从沿 x 轴的基线延伸的针状图。各个数据值由终止每个针状图的圆指示。
axis([0 data -2 2]);
grid on;
title('input binary code');
subplot(212)
%原始图像的频谱,fft快速傅里叶变换,abs绝对值和复数的模
plot(20*log(abs(fft(rand_data))));
axis([0 data -50 100]);
grid on;
title('spectrum of input binary code');
%% seriel to parallel %同时单极性码转为双极性码
for i=1:data
if rem(i,2)==1 %如果i是单数,rem表示取余数
if rand_data(i)==1%1为+1,0为-1
I(i)=1;%因为分为两路,所以此时一个码元占用原来两个码元的持续时间
I(i+1)=1;
else
I(i)=-1;
I(i+1)=-1;
end
else
if rand_data(i)==1
Q(i-1)=1;%因为i已经加到最大值了,只能累减
Q(i)=1;
else
Q(i-1)=-1;
Q(i)=-1;
end
end
end
figure(num_fig)%创建第二个窗口
num_fig = num_fig+1;
subplot(211);
%双极性码的频谱
plot(20*log(abs(fft(I))));
axis([0 data -50 140]);
grid on;
title('spectrum of I-channel data');
subplot(212);
plot(20*log(abs(fft(Q))));
axis([0 data -50 140]);
grid on;
title('spectrum of Q-channel data');
%% zero insertion 采样 ,此过程称为成形。成形的意思就是实现由消息到波形的转换,以便发射,脉冲成形应该是在基带调制之后。
zero=5; %sampling rate 25M HZ ,明白了,zero为过采样率,它等于 采样率fs/码速率。码元速率为 5MHZ。
%如果实际采样频率高于奈奎斯特频率,即为过采样
for i=1:zero*data % 采样点数目=过采样率*原码数目
if rem(i,zero)==1
Izero(i)=I(fix((i-1)/zero)+1); %fix,朝零方向四舍五入为最近的整数
Qzero(i)=Q(fix((i-1)/zero)+1);
else
Izero(i)=0;%插入零
Qzero(i)=0;
end
end
figure(num_fig)%第三个窗口
num_fig = num_fig+1;
subplot(211);
plot(20*log(abs(fft(Izero))));
axis([0 zero*data -20 140]);
grid on;
title('spectrum of I-channel after zero insertion');
subplot(212);
plot(20*log(abs(fft(Qzero))));
axis([0 zero*data -20 140]);
grid on;
title('spectrum of Q-channel after zero insertion');
%减小码间串扰
%pulse shape filter, 接着,将进行低通滤波,因为 随着传输速率的增大,基带脉冲的频谱将变宽
%如果不滤波(如升余弦滤波)进行低通滤波,后面加载频的时候可能会出现困难。
%% 平方根升余弦滤波器
% psf=rcosfir(rf,n_t,rate,fs,'sqrt') rate:过采样率,rf:滚降因子,n_t:滤波器阶数,fs:采样率
%用在调制或发送之前,用在解调或接受之后,用来降低过采样符号流带宽并不引发ISI(码间串扰)
%构造滤波器
NT=50;%滤波器阶数
N=2*zero*NT; % =500 zero为每个信号点采样的个数-1,插零个数,为什么这样算?
fs=25e6;%采样频率
rf=0.1;%滚降因子
psf=rcosfir(rf,NT,zero,fs,'sqrt');%psf大小为500
figure(num_fig)%第四张
num_fig = num_fig+1;
subplot(211);
plot(psf);
axis([200 300 -0.2 0.6]);
title('time domain response of pulse shaping filter');
grid on;
subplot(212);
plot(20*log(abs(fft(psf))));
axis([0 N -350 50]);
grid on;
title('transfer function of pulse shaping filter');
%用滤波器滤波
Ipulse=conv(Izero,psf);%conv 卷积
Qpulse=conv(Qzero,psf);
figure(num_fig)%第五张
num_fig = num_fig+1;
subplot(211);
plot(20*log(abs(fft(Ipulse))));
axis([0 zero*data+N -250 150]);%卷积后的长度zero*data+N-1
grid on;
title('spectrum of I-channel after impulse shaping filter');
subplot(212);
plot(20*log(abs(fft(Qpulse))));
axis([0 zero*data+N -250 150]);
grid on;
title('spectrum of Q-channel after impluse shaping filter');
%为什么数字信号传输也要过采样,成形滤波?
%答:过采样的数字信号处理起来对低通滤波器的要求相对较低,如果不过采样,滤波的时候滤波器需要很陡峭,指标会很严格
%成形滤波的作用是保证采样点不失真。如果没有它,那信号在经过带限信道后,眼图张不开,ISI非常严重。成形滤波的位置在基带调制之后。
%因为经成形滤波后,信号的信息已经有所损失,这也是为避免ISI付出的代价。换句话说,成形滤波的位置在载波调制之前,仅挨着载波调制。
%即:(发送端)插值(采样)-成形-滤波(LPF)-加载频(载波调制)-加噪声至(接收端)乘本振-低通-定时抽取-判决。
%% modulation调制
for i=1:zero*data+N %采样点数目改变 (因为卷积的 缘故)
t(i)=(i-1)/(fs); %这里因为假设载频与码速率大小相等,所以用载频fc乘以过采样率=采样率。时间要变化,1/fs是时间间隔
Imod(i)=Ipulse(i)*sqrt(2)*cos(2*pi*fc*t(i));%pi/4提出来了
Qmod(i)=Qpulse(i)*(-sqrt(2)*sin(2*pi*fc*t(i)));
end
sum=Imod+Qmod;
figure(num_fig)%6
num_fig = num_fig+1;
subplot(211);
plot(20*log(abs(fft(Imod))));
axis([0 zero*data+N -250 150]);
grid on ;
title('spectrum of I-channel after modulation');
subplot(212);
plot(20*log(abs(fft(Qmod))));
axis([0 zero*data+N -250 150]);
grid on;
title('spectrum of Q-channel after modulation');
figure(num_fig)%7
num_fig = num_fig+1;
plot(20*log(abs(fft(sum))));
axis([0 zero*data+N -250 150]);
grid on;
title('spectrum of sum after modulation');
figure(num_fig)%8
num_fig = num_fig+1;
plot(I,Q,'*'); %%只存在四种情况,故只显示四个点,分为了两路,一路是I路,一路是Q路,(0,1)、(1,0)、(1,1)、(0,0)
axis([-1.2 1.2 -1.2 1.2]);
grid on;
title('parallel constellation of sampler I and Q');
%% QPSK receiver
%demodulation解调
for i=1:zero*data+N
Idem(i)=sum(i)*sqrt(2)*cos(2*pi*fc*t(i));
Qdem(i)=sum(i)*(-sqrt(2)*sin(2*pi*fc*t(i)));
end
%matched filter匹配滤波器
mtf=rcosfir(rf,NT,zero,fs,'sqrt');%升余弦匹配滤波器
Imat=conv(Idem,mtf);
Qmat=conv(Qdem,mtf);
%data selection
for i=1:zero*data
Isel(i)=Imat(i+N);
Qsel(i)=Qmat(i+N);
end
%sampler %提取码元
for i=1:data
Isam(i)=Isel((i-1)*zero+1);
Qsam(i)=Qsel((i-1)*zero+1);
end
%decision threshold %判决门限
threshold=0.2;
for i=1:data
if Isam(i)>=threshold
Ifinal(i)=1;
else
Ifinal(i)=-1;
end
if Qsam(i)>=threshold
Qfinal(i)=1;
else
Qfinal(i)=-1;
end
end
%parallel to serial 并串转换
for i=1:data
if rem (i,2)==1
if Ifinal(i)==1
final(i)=1;
else
final(i)=0;
end
else
if Qfinal(i)==1
final(i)=1;
else
final(i)=0;
end
end
end
%% -----------------------------------------------------------------------------------------------------缁樺浘
figure(num_fig)%9
num_fig = num_fig+1;
subplot(221);
plot(20*log(abs(fft(Idem))));
axis([0 zero*data -200 150]);
grid on;
title('spectrum of I-channel after demodulation');
subplot(222);
plot(20*log(abs(fft(Qdem))));
axis([0 zero*data+N -200 150 ]);
grid on;
title('spectrum of Q-channel after demodulation');
subplot(223);
plot(20*log(abs(fft(Imat))));
axis([0 zero*data -400 200]);
grid on;
title('spectrum of I-channel after matched filter');
subplot(224);
plot(20*log(abs(fft(Qmat))));
axis([0 zero*data -400 200]);
grid on;
title('spectrum of Q-channel after matched filter');
figure(num_fig)%10
num_fig = num_fig+1;
subplot(221);
plot(20*log(abs(fft(Isam))));
axis([0 data -40 150]);
grid on;
title('spectrum of I-channel after sampler');
subplot(222);
plot(20*log(abs(fft(Qsam))));
axis([0 data -40 150 ]);
grid on;
title('spectrum of Q-channel after sampler');
subplot(223);
plot(20*log(abs(fft(Ifinal))));
axis([0 data -40 150]);
grid on;
title('spectrum of I-channel after decision threshold');
subplot(224);
plot(20*log(abs(fft(Qfinal))));
axis([0 data -40 150]);
grid on;
title('spectrum of Q-channel after decision threshold');
figure(num_fig)%11
num_fig = num_fig+1;
plot(Isel,Qsel);
axis([-1.6 1.6 -1.6 1.6]);
grid on;
title('constellation of matched filter output');%星座图
figure(num_fig)%12
num_fig = num_fig+1;
plot(Isam,Qsam,'X');
axis([-1.2 1.2 -1.2 1.2]);
grid on;
title('constellation of sampler');
figure(num_fig)%13
num_fig = num_fig+1;
subplot(211)
stem(final);
axis([0 data -2 2]);
grid on;
title('final received binary data');%获得的二进制数据
subplot(212)
plot(20*log(abs(fft(final))));
axis([0 data 0 100]);
grid on;
title('spectrum of final received binary data');
QPSK调制Matlab仿真
最新推荐文章于 2024-09-11 23:56:24 发布