QAM系统调制解调详解（带Matlab代码）

前言

早期，我整理了一个QPSK的调制解调过程，对其中的每一步进行公式推理和分析，能够说明一些问题，但是留了一个小坑说后续会提供matlab仿真一直没填。
近期，我需要对OFDM系统进行一些推导，我翻到的之前写的内容，顺便把过程又捋了一遍，编写了matlab代码（也算填坑了），希望能够更好的解释QAM调制解调过程，供大家一起交流学习。

QAM调制解调系统框图

QAM调制解调系统的公式推导

对于QAM系统，其实整体流程都是一致的，无非就是4QAM是每2个bit映射到一个复数，8QAM是3个bit，16QAM是4个比特，以此类推。。。
当然，随着阶数的不断升高，其抗干扰能力越弱，约容易出现误码，需要对系统进一步进行优化已保障调制解调的成果，本文仅为学习整体基本流程，不考虑系统优化等问题。

公式推导我们以16QAM为例，其他的大同小异。

1. 16QAM 4个bit用一个复数进行表示，我们先看一下16QAM的映射表和星座图。

QAM的调制

• Step1：假设我需要传4个bit 0001
• Step2：根据QAM的映射规则，对应复数是-3+j
• Step3：实部是 -3、虚部是 +1，ps虚部这儿图片上是错误的，应该是Img{}
• Step4：实部与本地载波cos相乘 -3cos(wt)、实部与本地载波sin相乘 +sin(wt)
• Step5：两部分相加即可得到发射信号
  重要的三角公式：$cos(a\pm b)=cosacosb\mp sinasinb$
  相加后计算结果为：$x(t)=3cos(wt)+sin(wt)=x(t)=-\surd 10cos(wt+arctan(1/3))$

无线信道

• Step6：无线信道我们暂时只考虑高斯白噪声，接收机将收到的信号可表示为$x(t)=-\surd 10cos(wt+arctan(1/3))+n$

QAM的解调

• Step7：接收信号分别与本振相乘，将分别得到
  重要的三角公式：
  $\cos A\cos B=\frac{1}{2}[\cos (A+B)+\cos (A-B)]$
  $\cos A\sin B=\frac{1}{2}[\sin (A+B)-\sin (A-B)]$
  $\sin A\sin B=\frac{1}{2}[\cos (A-B)-\cos (A+B)]$
  计算结果：
  $x(t)\ast cos(wt)=-\surd 10cos(wt+arctan(1/3))\ast cos(wt)=-3/2-3/2cos(2wt)+1/2sin(2wt)+noise$
  $x(t)\ast sin(wt)=-\surd 10cos(wt+arctan(1/3))\ast cos(wt)=1/2-1/2cos(2wt)-3/2sin(2wt)+noise$

• Step8：我们将Step7的公式进行拆分讲解，为啥要通过低通滤波器。
  x(t)∗cos(wt)==−3/2+noise−3/2cos(2wt)+1/2sin(2wt)
  上述结果中，我们看到乘以本振后，由三部分组成，即直流分量、高频分量cos、高频分量sin。我们需要并提供判断接收的依据是 直流分量 ，因此仅需要使用低通滤波器将高频分量滤除即可得到直流分量。

• Step9： 判决的方法有很多优化，我们只考虑最简单的判决，如下表

• Step10、11：判决完成后虚部先乘*j，再相加，得-3+j
• Step12：通过映射表可知-3+j，对应比特数据为0001

QAM调制解调系统的matlab详解

完整代码

第一部分，我先给出完整的matlab代码，没有画图部分，可以自行添加。

F0=100;
N=10;
Fs=F0*N;
Ts=1/Fs;
T0=1/F0;
n=(0:T0/Ts-1)*Ts;

%程序说明
%生成100组16QAM数据，并调制发送与接收

%step1 生成比特流 4*100 共400个0、1数据
bitdata = randi([0, 1], 1, 400);

%step2 映射到星座图，将每4个bit 转化为1个 a+jb形式
bitdata_IO = qammod(bitdata', 16 , InputType='bit',UnitAveragePower=false);

%step3.1 提取实部
bitdata_IO_real=real(bitdata_IO);

%step3.2 提取虚部
bitdata_IO_imag=imag(bitdata_IO);

%step4.1 实部与本地载波cos相乘
carrier_cos =cos(2*pi*F0*n);
bitdata_IO_real_cos = bitdata_IO_real * carrier_cos ;

%step4.2 虚部与本地载波sin相乘
carrier_sin =sin(2*pi*F0*n);
bitdata_IO_real_sin = bitdata_IO_imag * carrier_sin ;

%step5 实部与虚部相加发射
TR = (bitdata_IO_real_cos + bitdata_IO_real_sin)' ;

%step6 通过信道，增加噪声
RX = awgn(TR,10);


%step7.1 接收信号与本地载波cos相乘
RX_cos = RX.*carrier_cos';

%step7.2 接收信号与本地载波sin相乘
RX_sin = RX.*carrier_sin';


%step8.1 实部通过低通滤波器
Fc = 10; % 通带截止频率 50 Hz
Wn = Fc / (Fs / 2);
hh = fir1(10, Wn, hann(10 + 1));% 使用汉宁窗设计 FIR 低通滤波器
y_RX_cos = filter(hh, 1, RX_cos);

%step8.2 虚部通过低通滤波器
y_RX_sin = filter(hh, 1, RX_sin);

%step9.1 实部抽样判决通过低通滤波器
y_RX_cos_mean = mean(y_RX_cos);
y_RX_cos_dec = y_RX_cos_mean;
y_RX_cos_dec(y_RX_cos_dec>max(y_RX_cos_mean)/2)=3;
y_RX_cos_dec(y_RX_cos_dec>0 & y_RX_cos_dec<max(y_RX_cos_mean)/2)=1;
y_RX_cos_dec(y_RX_cos_dec<min(y_RX_cos_mean)/2)=-3;
y_RX_cos_dec(y_RX_cos_dec<0 & y_RX_cos_dec>min(y_RX_cos_mean)/2)=-1;

%step9.2 虚部抽样判决通过低通滤波器
y_RX_sin_mean = mean(y_RX_sin);
y_RX_sin_dec = y_RX_sin_mean;
y_RX_sin_dec(y_RX_sin_dec>max(y_RX_sin_mean)/2)=3;
y_RX_sin_dec(y_RX_sin_dec>0 & y_RX_sin_dec<max(y_RX_sin_mean)/2)=1;
y_RX_sin_dec(y_RX_sin_dec<min(y_RX_sin_mean)/2)=-3;
y_RX_sin_dec(y_RX_sin_dec<0 & y_RX_sin_dec>min(y_RX_sin_mean)/2)=-1;

%step10、11 虚部与实部相加，还原复数
y_and = y_RX_cos_dec + 1i.*y_RX_sin_dec;

%step12 虚部与实部相加，还原复数
y = qamdemod(y_and, 16,OutputType='bit',UnitAveragePower=false);
y = reshape(y,1,400);

[number,ratio] = biterr(y,bitdata);
disp(['误码数：', num2str(number)]);
disp(['误码率：', num2str(ratio)]);

代码详解

step1：生成比特流

bitdata = randi([0, 1], 1, 400);

说明：代码中生成了400个bit，经过16QAM调制，每4位映射1个复数，共100组。

step2：映射到星座图

bitdata_IO = qammod(bitdata', 16 , InputType='bit',UnitAveragePower=false);
scatterplot(bitdata_IO);

说明：按照映射表进行映射。

step3：映射到星座图

%step3.1 提取实部
bitdata_IO_real=real(bitdata_IO);

%step3.2 提取虚部
bitdata_IO_imag=imag(bitdata_IO);

说明：没啥好说的，取就行了

step4：与本地载波相乘

%step4.1 实部与本地载波cos相乘
carrier_cos =cos(2*pi*F0*n);
bitdata_IO_real_cos = bitdata_IO_real * carrier_cos ;

%step4.2 虚部与本地载波sin相乘
carrier_sin =sin(2*pi*F0*n);
bitdata_IO_real_sin = bitdata_IO_imag * carrier_sin ;

figure(2)
subplot(2,1,1)
plot(n,bitdata_IO_real_cos)
title('实部与本地载波cos相乘');
xlabel('幅度');
ylabel('时间');
grid on;
subplot(2,1,2)
plot(n,bitdata_IO_real_sin)
title('虚部与本地载波sin相乘');
xlabel('幅度');
ylabel('时间');
grid on;
subplot(2,2,3);
fft_bitdata_IO_real_cos = fft(bitdata_IO_real_cos',10*length(bitdata_IO_real_cos));
fft_bitdata_IO_real_cos_shift = fftshift(fft_bitdata_IO_real_cos);
f=(-length(fft_bitdata_IO_real_cos_shift)/2:length(fft_bitdata_IO_real_cos_shift)/2-1).*(Fs/length(fft_bitdata_IO_real_cos_shift));
plot(f,abs(fft_bitdata_IO_real_cos_shift));
title('频域：实部与本地载波cos相乘');
xlabel('幅度');
ylabel('频率');
grid on;
subplot(2,2,4);
fft_bitdata_IO_real_sin = fft(bitdata_IO_real_sin',10*length(bitdata_IO_real_sin));
fft_bitdata_IO_real_sin_shift = fftshift(fft_bitdata_IO_real_sin);
f=(-length(fft_bitdata_IO_real_sin_shift)/2:length(fft_bitdata_IO_real_sin_shift)/2-1).*(Fs/length(fft_bitdata_IO_real_sin_shift));
plot(f,abs(fft_bitdata_IO_real_sin_shift));
title('频域：虚部与本地载波sin相乘');
xlabel('幅度');
ylabel('频率');
grid on;

说明：
可以看到，实部与虚部分别用cos与sin相乘，出现了四种波形。
同时，与本地载波相乘，也都般频到了载波上（100Hz）
tip：由于频域计算时取得是abs绝对值，因此频谱上只能看到2种，其实是正负4种。

step5：实部与虚部相加发射

%step5 实部与虚部相加发射
TR = (bitdata_IO_real_cos + bitdata_IO_real_sin)' ;
figure(3)
subplot(2,1,1)
plot(n,TR)
title('时域：实部与虚部相加发射');
xlabel('幅度');path
ylabel('时间');
grid on;
subplot(2,1,2);
fft_TR = fft(TR,10*length(TR));
fft_TR_shift = fftshift(fft_TR);
f=(-length(fft_TR_shift)/2:length(fft_TR_shift)/2-1).*(Fs/length(fft_TR_shift));
plot(f,abs(fft_TR_shift));
title('频域：实部与虚部相加发射');
xlabel('幅度');
ylabel('频率');
grid on;

说明：经过实部与虚部得组合，一共可得到8个幅值和相位都不相同得波形，从频谱上也可以看出。

step6：通过AWGN信道

RX = awgn(TR,10);

figure(3)
subplot(2,1,1)
plot(n,RX)
title('时域：通过AWGN信道，增加噪声');
xlabel('幅度');path
ylabel('时间');
grid on;
subplot(2,1,2);
fft_RX = fft(RX,10*length(RX));
fft_RX_shift = fftshift(fft_RX);
f=(-length(fft_RX_shift)/2:length(fft_RX_shift)/2-1).*(Fs/length(fft_RX_shift));
plot(f,abs(fft_RX_shift));
title('频域：通过AWGN信道，增加噪声');
xlabel('幅度');
ylabel('频率');
grid on;

说明：增加噪声后，100组波形出现了不同，我添加的是10dB噪声，也可以尝试其他的。

step7：接收信号与本振相乘

%step7.1 接收信号与本地载波cos相乘
RX_cos = RX.*carrier_cos';
%step7.2 接收信号与本地载波sin相乘
RX_sin = RX.*carrier_sin';


figure(5)
subplot(2,2,1)
plot(n,RX_cos)
title('时域：接收信号与本振cos相乘');
xlabel('幅度');path
ylabel('时间');
grid on;
subplot(2,2,2)
plot(n,RX_sin)
title('时域：接收信号与本振sin相乘');
xlabel('幅度');
ylabel('时间');
grid on;
subplot(2,2,3);
fft_RX_cos = fft(RX_cos,10*length(RX_cos));
fft_RX_cos_shift = fftshift(fft_RX_cos);
f=(-length(fft_RX_cos_shift)/2:length(fft_RX_cos_shift)/2-1).*(Fs/length(fft_RX_cos_shift));
plot(f,abs(fft_RX_cos_shift));
title('频域：接收信号与本振cos相乘');
xlabel('幅度');
ylabel('频率');
grid on;
subplot(2,2,4);
fft_RX_sin = fft(RX_sin,10*length(RX_sin));
fft_RX_sin_shift = fftshift(fft_RX_sin);
f=(-length(fft_RX_sin_shift)/2:length(fft_RX_sin_shift)/2-1).*(Fs/length(fft_RX_sin_shift));
plot(f,abs(fft_RX_sin_shift));
title('频域：接收信号与本振sin相乘');
xlabel('幅度');
ylabel('频率');
grid on;

说明：通过与本振相乘，从频谱上可以考到，每个信号都出现直流分量（0Hz）和高频分量（200Hz）

step8：通过低通滤波器，滤出高频分量

%step8.1 实部通过低通滤波器
Fc = 10; % 通带截止频率 50 Hz
Wn = Fc / (Fs / 2);
hh = fir1(10, Wn, hann(10 + 1));% 使用汉宁窗设计 FIR 低通滤波器
freqz(hh,1,512,Fs)
y_RX_cos = filter(hh, 1, RX_cos);
%step8.1 虚部通过低通滤波器
y_RX_sin = filter(hh, 1, RX_sin);

figure(6)
subplot(2,2,1)
plot(n,y_RX_cos)
title('时域：与cos相乘后通过低通滤波器');
xlabel('幅度');path
ylabel('时间');
grid on;
subplot(2,2,2)
plot(n,y_RX_sin)
title('时域：与sin相乘后通过低通滤波器');
xlabel('幅度');
ylabel('时间');
grid on;
subplot(2,2,3);
fft_y_RX_cos = fft(y_RX_cos,10*length(y_RX_cos));
fft_y_RX_cos_shift = fftshift(fft_y_RX_cos);
f=(-length(fft_y_RX_cos_shift)/2:length(fft_y_RX_cos_shift)/2-1).*(Fs/length(fft_y_RX_cos_shift));
plot(f,abs(fft_y_RX_cos_shift));
title('频域：与sin相乘后通过低通滤波器');
xlabel('幅度');
ylabel('频率');
grid on;
subplot(2,2,4);
fft_y_RX_sin = fft(y_RX_sin,10*length(y_RX_sin));
fft_y_RX_sin_shift = fftshift(fft_y_RX_sin);
f=(-length(fft_y_RX_sin_shift)/2:length(fft_y_RX_sin_shift)/2-1).*(Fs/length(fft_y_RX_sin_shift));
plot(f,abs(fft_y_RX_sin_shift));
title('频域：接收信号与本振sin相乘');
xlabel('幅度');
ylabel('频率');
grid on;

说明：
滤波器实在业余，随便弄了一个滤波器，截至频率是50Hz，滤波器特性可以看响应图。这部分应该是可以优化的，缓冲带越窄越好。
通过滤波器后，我们可以看到基本上只保留的直流部分，其他频率的信号都被抑制了。

step9：抽样判决

%step9.1 实部抽样判决通过低通滤波器
y_RX_cos_mean = mean(y_RX_cos);
y_RX_cos_dec = y_RX_cos_mean;
y_RX_cos_dec(y_RX_cos_dec>max(y_RX_cos_mean)/2)=3;
y_RX_cos_dec(y_RX_cos_dec>0 & y_RX_cos_dec<max(y_RX_cos_mean)/2)=1;
y_RX_cos_dec(y_RX_cos_dec<min(y_RX_cos_mean)/2)=-3;
y_RX_cos_dec(y_RX_cos_dec<0 & y_RX_cos_dec>min(y_RX_cos_mean)/2)=-1;

%step9.2 虚部抽样判决通过低通滤波器
y_RX_sin_mean = mean(y_RX_sin);
y_RX_sin_dec = y_RX_sin_mean;
y_RX_sin_dec(y_RX_sin_dec>max(y_RX_sin_mean)/2)=3;
y_RX_sin_dec(y_RX_sin_dec>0 & y_RX_sin_dec<max(y_RX_sin_mean)/2)=1;
y_RX_sin_dec(y_RX_sin_dec<min(y_RX_sin_mean)/2)=-3;
y_RX_sin_dec(y_RX_sin_dec<0 & y_RX_sin_dec>min(y_RX_sin_mean)/2)=-1;

说明：通过第八步得频谱可以看到，虽然已经滤除了大部分非直流量，但是还存在部分其他频率。为了能得到直流分量，可以求均值。

step10、11：虚部与实部相加，还原复数

%step10 虚部与实部相加，还原复数
y_and = y_RX_cos_dec + 1i.*y_RX_sin_dec;

说明：没啥好说的，正常合并

step12：复数按照映射，还原为bit

y = qamdemod(y_and, 16,OutputType='bit',UnitAveragePower=false);
y = reshape(y,1,400);

[number,ratio] = biterr(y,bitdata);
disp(['误码数：', num2str(number)]);
disp(['误码率：', num2str(ratio)]);

说明：我们可以与发送的数据流进行对比，查看误码数和误码率。当然，如果信噪比过小，必定会出现误码
输出结果：
误码数：0
误码率：0

搞定，填坑结束，接下来可能写一下OFDM的文章，但愿有动力！！！！！！！！