MATLAB仿真:FSK调制解调

版权声明:本文为博主原创文章,未经博主允许不得转载。

本文档详细介绍如何使用MATLAB仿真一个简单但完整的通信系统——FSK调制解调。

FSK调制解调框图如下图所示:

fskmod_demod_diagram

系统分析:

  • 该系统的输入信号是基带信号。在此处,一定要清楚该基带信号是数字信号,且该数字信号是对模拟信号数字化(ADC)得来的。模数转换的流程是:采样->量化->编码。在采样的过程中,采样率必须满足采样定理。经过模数转换之后,得到了一串基带信号,就是该系统的输入信号。
  • 载波频率为f1和f2
  • 信道为加性高斯白噪声(Additive White Gauss Noise, AWGN)信道
  • 带通滤波器BPF1负责获得载频为f1的调制信号,带通滤波器BPF2负责获得载频为f2的调制信号
  • 相乘器的目的是实现相干解调。相乘器的输入信号只有一根谱线f1(或f2),输出有两根谱线0和2f1(或0和2f2)
  • 低通滤波器LPF只允许零频附近信号通过,高频信号2f1和2f2无法通过。通过滤波器之后就可以得到基带信号
  • 抽样判决器对LPF输出的基带信号进行判决

Note: 在模数转换过程中,一个采样点经过量化编码之后可能会被编码成多个比特。

 

系统配置

  • 采样率:2000Hz(Note: 此处的采样率与上文讲述的模数转换过程中的采样率不一样,这里的采样率是为了让画出来的图看起来是连续的。比如说,有一串基带信号为signal = [0 1 0 1 1 0 1],如果用 plot(signal)画基带信号波形,那么该波形是惨不忍睹的,如果每个01都画2000个,那么,基带波形就是比较符合预期的。
  • 模拟发送的比特数:20
  • 每个bit持续的时间:1s
  • 载频1的模拟频率:60Hz
  • 载频2的模拟频率:300Hz
  • 信噪比:2dB
  • 数字滤波器通带最大波动:0.1dB
  • 数字滤波器阻带最小衰减:60dB
  • 数字滤波器过渡带宽:0.01pi

有关数字滤波器设置的有关问题,请移步:Digital Signal Process – 数字滤波器的设计

下面只有代码,想看代码和仿真结果的请访问:http://www.wydfx.xyz/wp-dsp/fsk_mod_demod.html

%--------------------------------------------------------------------------
% Function    : FSK_main.m
% Description : A MATLAB program to simulate FSK modulation and
%               demoudulation.
% Input       : None
% Output      : None
% Version     : 1.05
% Author      : chengkai
% History     : When         Who                 What
% structure
%                                                
%--------------------------------------------------------------------------
clc,clear
close all
% 系统配置
sampleRate = 2000;                          % 采样率
dt = 1 / sampleRate;                        % 采样时间间隔
SimBitsNum = 20;                            % 模拟比特数
bitTime = 1;                                % 1bit持续的时间
bitSampleNum = sampleRate * bitTime;        % 1bit的采样点数
totalTime = bitTime * SimBitsNum;           % 总时间
totalSampleNum = sampleRate * totalTime;    % 采样点的个数
t = 0:dt:(totalTime - dt);                  % SimBitsNum个比特所持续的时间
t = t';                                     % 转置成列向量'
f1 = 60;                                    % 载频1的模拟频率(Hz)
f2 = 300;                                   % 载频2的模拟频率(Hz)
w1 = 2 * pi * f1 * dt / pi;                 % 载频1的数字频率/pi 在设计FIR滤波器时使用
w2 = 2 * pi * f2 * dt / pi;                 % 载频2的数字频率/pi 在设计FIR滤波器时使用
SNR = 2;                                    % 信噪比(单位:dB)
ap = 0.1;                                   % 通带最大波动
as = 60;                                    % 阻带最小衰减
tr_width = 0.01 * pi;                       % 数字滤波器过渡带宽
% 产生发送比特
% sendBits = randi([0 1],[1, SimBitsNum]);
sendBits = round(rand(1, SimBitsNum));              % 模拟发送的比特
sendBitsReverse = ~sendBits;
% 抽样
% 注意:此处使用了矩阵相乘
sendBitsSample = (ones(1, bitSampleNum))' * sendBits; %'
sendBitsSample = sendBitsSample(:);                 % 转换成列向量
sendBitsReverseSample = (ones(1, bitSampleNum))' * sendBitsReverse; %'
sendBitsReverseSample = sendBitsReverseSample(:);   % 转换成列向量

% FSK调制
fsk1 = sendBitsSample .* cos(2 * pi * f1 .* t);
fsk2 = sendBitsReverseSample .* cos(2 * pi * f2 .* t);
fsk = fsk1 + fsk2;
subplot(211)
plot(fsk)
axis([0 totalSampleNum min(fsk)-0.2 max(fsk)+0.2]);
grid on
title('FSK调制之后的时域波形')

% 添加噪声
% ====== AWGN noise ====== 
noisePower  = 10^(-SNR/10);
noise       = sqrt(noisePower) * (randn(length(fsk), 1));
fskSigWithNoise = fsk + noise;
subplot(212)
plot(fskSigWithNoise)
axis([0 totalSampleNum min(fskSigWithNoise)-0.2 max(fskSigWithNoise)+0.2]);
grid on
title('添加噪声之后的时域波形')

% 绘制添加噪声前后的信号幅度谱
figure
subplot(211)
plot((0:totalSampleNum-1)*(2/totalSampleNum), abs(fft(fsk)) / sqrt(totalSampleNum))
grid on
title('未加噪声的FSK调制信号的幅度谱')
subplot(212)
plot((0:totalSampleNum-1)*(2/totalSampleNum), abs(fft(fskSigWithNoise)) / sqrt(totalSampleNum))
grid on
title('添加噪声的FSK调制信号的幅度谱')

% 设计两个带通数字滤波器
hammingWinN = ceil(6 * pi / tr_width) + 1;
b1 = fir1(hammingWinN, [max(w1-0.02, 0) w1+0.02]);
b2 = fir1(hammingWinN, [w2-0.02 w2+0.02]);
freqz( dfilt.dffir(b1));
title('f1滤波器')
freqz( dfilt.dffir(b2));
title('f2滤波器')

% 滤波
H1 = filter(b1, 1, fskSigWithNoise);
H2 = filter(b2, 1, fskSigWithNoise);
figure
subplot(211)
plot(t, H1);
grid on
xlabel('t')
ylabel('幅度')
title('滤波之后的时域波形(fc = 20Hz)')
subplot(212)
plot(t, H2);
grid on
xlabel('t')
ylabel('幅度')
title('滤波之后的时域波形(fc = 120Hz)')

% 解调
% 使用相干解调的方法
sw1 = H1 .* H1;
sw2 = H2 .* H2;

figure
subplot(211)
plot(t, sw1)
title('经过相乘器H1后的波形')
xlabel('t')
ylabel('幅度')
grid on

subplot(212)
plot(t, sw2)
title('经过相乘器H2后的波形')
xlabel('t')
ylabel('幅度')
grid on

% 相乘之后的幅度谱
figure
subplot(211)
plot((0:totalSampleNum-1)*(2/totalSampleNum), abs(fft(sw1)) / sqrt(totalSampleNum))
grid on
title('sw1x')
subplot(212)
plot((0:totalSampleNum-1)*(2/totalSampleNum), abs(fft(sw2)) / sqrt(totalSampleNum))
grid on
title('sw2x')

% 设计低通滤波器
hammingWinN = ceil(6 * pi / (0.001 * pi)) + 1;
% LPF = fir1(hammingWinN, 0.05);
% freqz( dfilt.dffir(LPF));
% title('低通滤波器')
LPF = fir1(hammingWinN, 0.005);
% LPF = fir1(101, [2/800 10/800]);
freqz(dfilt.dffir(LPF));
% 经过低通滤波器
title('低通滤波器')
st1 = filter(LPF, 1, sw1);
st2 = filter(LPF, 1, sw2);
figure
subplot(211)
plot(t, st1)
title('经过低通滤波器后的波形1')
xlabel('t')
ylabel('幅度')
grid on

subplot(212)
plot(t, st2)
title('经过低通滤波器后的波形1')
xlabel('t')
ylabel('幅度')
grid on

% 判决
delay = 1.5;        % 由滤波器产生的时延,单位为符号,该数值是从图上直接观察到的
% 确定判决时刻
judgePoint = (bitSampleNum/2 + delay * sampleRate : bitSampleNum : totalSampleNum);
judgeFinal = (st1(judgePoint) > 0.3)'; %'
judgeNum = length(judgeFinal);      % 最后可以完整判断的符号的数量
errorRate = length(find(xor(sendBits(1 : judgeNum), judgeFinal) == 1)) / judgeNum;
fprintf('系统误码率为 %d\n', errorRate)