【信号识别】基于matlab深度学习CNN信号调制分类【含Matlab源码 2066期】

⛄一、深度学习CNN信号调制分类概述

1 背景介绍
在通信信号处理领域, 特别是在非协作通信信号盲解调研究领域, 每时隙突发信号的调制方式不同, 必须进行信号的调制方式自动识别。信号的调制方式识别效果会影响整个盲解调系统。传统的专家系统已经被研究多年, 然而依然存在准确性低、效率低等问题。近年来, 深度学习方法广泛应用于图像处理[1,2]和语音识别[3,4]领域, 取得了不错的效果。信号星座图是表征数字调制信号的一种重要方式, 它能够直观地将信号以图片几何特征的形式展示出来, 便于人们观察和分析。传统的人眼分辨方法依赖于一定的专业知识、耗时且识别准确率较低。我们将深度学习方法引用到通信信号处理领域, 提出了一种基于卷积神经网络的信号调制方式自动识别方法。

2 生成数据集
2.1原始数据的生成
实验需要获取多种调制方式、宽信噪比范围、带有标签的数据, 这在短时间内难以通过采集真实信号来实现。本次实验所需数据通过仿真产生, 实验选取部分数据进行训练学习, 使用另外的数据进行验证测试。仿真生成数据过程主要有生成随机0-M数字 (M为调制阶数) 、MQAM/MPSK调制、成形滤波、载波上变频、信道仿真 (加性高斯白噪声信道) 、IQ正交解调、匹配滤波、重构信号星座图并保存。整个生成数据流程如下所示。

图1 生成信号流程图

本次仿真共生成7种调制方式 (BPSK、QPSK、8PSK、QAM8、QAM16、QAM32、QAM64) 信号, 每种调制方式信号信噪比为-20dB、-18dB、…、+16dB、+18dB。每种调制方式在某一信噪比下共有1000条数据, 生成1000张星座图片, 每张图片包含1000符号信息 (即每张星座图图上有1000个点) 。仿真生成的基带信号经成形滤波后时域图形如下所示。

图2 经成形滤波后的信号时域图

2.2 IQ正交解调和星座图
在通信系统接收方通常需要将接收到的中频信号变成基带信号这就需要数字下变频操作。数字下变频 (DDC) 通常包括IQ正交解调和低通滤波两部分。IQ正交解调需要知道信号的载波频率, 低通滤波器截止频率的设置和信号带宽相关。经过低通滤波过后即可得到IQ两路信号, 此时得到的IQ路信号还需要经过匹配滤波、下采样符号同步等过程, 才能重构信号星座图。

图3 数字下变频原理说明图

图4 IQ两路信号和幅度相位关系说明图

IQ信号能够同时包含信号的幅度和相位等信息, IQ信号可以轻松地形成一个复信号。大多数数字调制将数据映射到IQ平面中的许多离散点, 称之为星座图。

3 基于CNN的调制方式识别工作原理
3.1卷积网络结构
基于CNN的调制方式识别模型网络框图如下所示。输入层输入的是150×120像素尺寸大小的信号星座图, 经过3个卷积层、3个池化层, 扁平层, 到最后按照7种调制类型输出。训练集数据都是带有标签的, 即事先知道某信号的调制方式, 不断学习图片对应调制类型, 整个系统通过梯度反向传播机制, 不断调节其中的相关参数, 使整个系统的输出结果和标签结果的误差最小。这就是整个模型系统的一个学习过程。

图5 CNN网络模型结构图

图6 卷积操作池化说明图
卷积池化操作是CNN中重要的组成部分, 上图左图中假设了4×5大小的输入矩阵数据, 和3×3的卷积核进行卷积运算的过程, 得到信号的特征图。上图右图是4×4的特征层进行2×2最大池化的运算过程。

CNN的功能是从特定模型中提取特征, 然后根据特征进行分类识别、预测或做出决策。在此过程中, 最重要的一步是特征提取, 即如何提取能够最好地区分事物的特征。如果提取的特征不能进行分类, 那么特征提取步骤将毫无意义。网络模型中卷积层层数越多, 更容易把握输入信号的细微特征[6]。然而, 在深度神经网络的设计中, 应该在每个阶段考虑卷积层数和核大小, 尝试以最少的计算量获得最佳结果, 网络设计需要平衡网络结构的宽度和深度。

对于相同的CNN网络结构, 迭代的次数, 数据量的大小和学习率都会影响模型的分类结果和有效性。这些参数的设置都需要在实践中不断尝试、不断分析修改。在本实验中, 设计的CNN模型使用3层卷积网络层, 卷积核大小为3×3, 迭代次数为40, 随着训练量的变化调整参数大小。

3.2每层输入和输出大小
输入层:150×120像素大小的信号星座图。卷积层1 (Conv1) :由3×3内核生成的148×118个特征图。池化层第1层 (Pooling1) :从2×2区域进行二次采样后的74×59个像素图。卷积层2 (Conv2) :由3×3内核生成的72×57个特征映射。池化层第2层 (Pooling2) :从2×2区域进行二次采样后的36×28个特征映射。卷积层3 (Conv3) :由3×3内核生成的34×26个特征映射。池化第3层 (Pooling3) :从2×2区域进行二次采样后的17×13个特征映射。完全连接层 (F1atten) :从Pooling 3的所有像素转换的14144个节点, 即17×13×64=14144。Dense1:64, Dense2 (输出层Output) 7个结点, Softmax分类器分类输出, 输出7个节点, 代表七种不同的调制方法。

3.3星座图
由于卷积和池化操作, CNN可以从图像细微差别中检测到信号特征, 来分类识别信号。我们使用七种最常用的调制方法 (BPSK, QPSK, 8PSK, QAM8, QAM16, QAM32, QAM64) 作为例子来说明。从下图可以看出, BPSK信号星座图分布在x=0轴两端, 成2个块团状;QPSK信号分布X轴Y轴正负两侧。呈现4个团块状;8PSK分布成8个块团, 且每个块团离坐标原点的距离相等, 8个块团形成了一个圆。QAM16均匀分布在坐标轴四个象限, 每个现象有四个团块, 共16个团块。综上所述, 不同数字调制信号星座图的几何形状不同, 在实际工程中, 通信工程师经常会根据星座图形状来判断信号的调制方式。

图7 在SNR为10dB时BPSK、QPAK、8PSK、QAM8、QAM16、QAM64调制信号的星座图

⛄二、部分源代码

clear all;
clc;

%% 1 生成用于训练的波形
% 为每种调制类型生成 10000 个帧，其中 80% 用于训练，10% 用于验证，10% 用于测试。我们在网络训练阶段使用训练和验证帧。使用测试帧获得最终分类准确度。每帧的长度为 1024 个样本，采样率为 200 kHz。对于数字调制类型，八个样本表示一个符号。网络根据单个帧而不是多个连续帧（如视频）作出每个决定。假设数字和模拟调制类型的中心频率分别为 900 MHz 和 100 MHz。
% 要快速运行此示例，请使用经过训练的网络并生成少量训练帧。要在计算机上训练网络，请选择 “train network now” 选项(即将 trainNow 设置为true)。
trainNow = false;
if trainNow == true
numFramesPerModType = 100;
else
numFramesPerModType = 100;
end
percentTrainingSamples = 80;
percentValidationSamples = 10;
percentTestSamples = 10;

sps = 8; % Samples per symbol
spf = 1024; % Samples per frame
symbolsPerFrame = spf / sps;
fs = 200e3; % Sample rate
fc = [902e6 100e6]; % Center frequencies

%% 1.1 信道创建
% 让每帧通过信道并具有
% AWGN
% 莱斯多径衰落
% 时钟偏移，导致中心频率偏移和采样时间偏移
% 由于本示例中的网络基于单个帧作出决定，因此每个帧必须通过独立的信道

%% 1.1.1 AWGN
% 添加 SNR 为 30 dB 的 AWGN。由于帧经过归一化，因此噪声标准差可以计算为
SNR = 30;
std = sqrt(10.^(-SNR/10))
% 使用 comm.AWGNChannel 实现
awgnChannel = comm.AWGNChannel(…
‘NoiseMethod’, ‘Signal to noise ratio (SNR)’, …
‘SignalPower’, 1, …
‘SNR’, SNR);

%% 1.1.2 莱斯多径衰落
% 使用 comm.RicianChannel System object 实现通过莱斯多径衰落信道。假设延迟分布为 [0 1.8 3.4] 个样本，对应的平均路径增益为 [0 -2 -10] dB。K 因子为 4，最大多普勒频移为 4 Hz，等效于 900 MHz 的步行速度。使用以下设置实现信道。
multipathChannel = comm.RicianChannel(…
‘SampleRate’, fs, …
‘PathDelays’, [0 1.8 3.4]/fs, …
‘AveragePathGains’, [0 -2 -10], …
‘KFactor’, 4, …
‘MaximumDopplerShift’, 4);

%% 1.1.3 时钟偏移
% .时钟偏移是发送器和接收器的内部时钟源不准确造成的。时钟偏移导致中心频率（用于将信号下变频至基带）和数模转换器采样率不同于理想值。信道仿真器使用时钟偏移因子 C，表示为 C=1+Δclock106，其中 Δclock 是时钟偏移。对于每个帧，通道基于 [?maxΔclock maxΔclock] 范围内一组均匀分布的值生成一个随机 Δclock 值，其中 maxΔclock 是最大时钟偏移。时钟偏移以百万分率 (ppm) 为单位测量。对于本示例，假设最大时钟偏移为 5 ppm。
maxDeltaOff = 5;
deltaOff = (rand()2maxDeltaOff) - maxDeltaOff;
C = 1 + (deltaOff/1e6);
% 频率偏移
% 基于时钟偏移因子 C 和中心频率，对每帧进行频率偏移。使用 comm.PhaseFrequencyOffset 实现信道。
offset = -(C-1)*fc(1);
frequencyShifter = comm.PhaseFrequencyOffset(…
‘SampleRate’, fs, …
‘FrequencyOffset’, offset);
% 采样率偏移
% 基于时钟偏移因子 C，对每帧进行采样率偏移。使用 interp1 函数实现通道，以 C×fs 的新速率对帧进行重新采样。

%% 1.1.4 合成后的信道
% 使用 helperModClassTestChannel 对象对帧应用所有三种信道衰落
channel = helperModClassTestChannel(…
‘SampleRate’, fs, …
‘SNR’, SNR, …
‘PathDelays’, [0 1.8 3.4] / fs, …
‘AveragePathGains’, [0 -2 -10], …
‘KFactor’, 4, …
‘MaximumDopplerShift’, 4, …
‘MaximumClockOffset’, 5, …
‘CenterFrequency’, 902e6)
% 您可以使用 info 对象函数查看有关通道的基本信息。
chInfo = info(channel);

%% 1.2 波形生成
% 创建一个循环，它为每种调制类型生成信道衰落的帧，并将这些帧及其对应标签存储在 frameStore 中。从每帧的开头删除随机数量的样本，以去除瞬变并确保帧相对于符号边界具有随机起点。
% Set the random number generator to a known state to be able to regenerate
% the same frames every time the simulation is run

tic

modulationTypes = categorical([“BPSK”, “QPSK”, “8PSK”, …
“16QAM”, “64QAM”, “PAM4”, “GFSK”, “CPFSK”, …
“B-FM”, “DSB-AM”, “SSB-AM”]);

numModulationTypes = length(modulationTypes);

channelInfo = info(channel);
frameStore = helperModClassFrameStore(…
numFramesPerModType*numModulationTypes,spf,modulationTypes);
transDelay = 50;

for modType = 1:numModulationTypes
fprintf(‘%s - Generating %s frames\n’, …
datestr(toc/86400,‘HH:MM:SS’), modulationTypes(modType))
numSymbols = (numFramesPerModType / sps);
dataSrc = getSource(modulationTypes(modType), sps, 2*spf, fs);
modulator = getModulator(modulationTypes(modType), sps, fs);
if contains(char(modulationTypes(modType)), {‘B-FM’,‘DSB-AM’,‘SSB-AM’})
% Analog modulation types use a center frequency of 100 MHz
channel.CenterFrequency = 100e6;
else
% Digital modulation types use a center frequency of 902 MHz
channel.CenterFrequency = 902e6;
end

for p=1:numFramesPerModType
% Generate random data
x = dataSrc();

% Modulate
y = modulator(x);% Pass through independent channels
rxSamples = channel(y);% Remove transients from the beginning, trim to size, and normalize
frame = helperModClassFrameGenerator(rxSamples, spf, spf, transDelay, sps);% Add to frame store
add(frameStore, frame, modulationTypes(modType));

end
end

% 接下来，将帧分为训练数据、验证数据和测试数据。默认情况下，frameStore 将 I/Q 基带样本按行放置在输出帧中。输出帧的大小为 [2xspf×1×N]，其中第一行是同相采样，第二行是正交采样。
[mcfsTraining,mcfsValidation,mcfsTest] = splitData(frameStore,…
[percentTrainingSamples,percentValidationSamples,percentTestSamples]);
[rxTraining,rxTrainingLabel] = get(mcfsTraining);
[rxValidation,rxValidationLabel] = get(mcfsValidation);
[rxTest,rxTestLabel] = get(mcfsTest);

size(rxTraining)

% Plot the amplitude of the real and imaginary parts of the example frames
% against the sample number
plotTimeDomain(rxTest,rxTestLabel,modulationTypes,fs)

% Plot a spectrogram of the example frames
plotSpectrogram(rxTest,rxTestLabel,modulationTypes,fs,sps)

% 通过确保标签（调制类型）分布均匀，避免训练数据中的类不平衡。绘制标签分布图，以检查生成的标签是否分布均匀。
% Plot the label distributions
figure
subplot(3,1,1)
histogram(rxTrainingLabel)
title(“Training Label Distribution”)
subplot(3,1,2)
histogram(rxValidationLabel)
title(“Validation Label Distribution”)
subplot(3,1,3)
histogram(rxTestLabel)
title(“Test Label Distribution”)

%% 2 训练 CNN
% 本示例使用的 CNN 由六个卷积层和一个全连接层组成。除最后一个卷积层外，每个卷积层后面都有一个批量归一化层、修正线性单元 (ReLU) 激活层和最大池化层。在最后一个卷积层中，最大池化层被一个平均池化层取代。输出层具有 softmax 激活。有关网络设计指导原则，请参阅Deep Learning Tips and Tricks。
dropoutRate = 0.5;
numModTypes = numel(modulationTypes);
netWidth = 1;
filterSize = [1 sps];
poolSize = [1 2];
modClassNet = [
imageInputLayer([2 spf 1], ‘Normalization’, ‘none’, ‘Name’, ‘Input Layer’)

convolution2dLayer(filterSize, 16*netWidth, ‘Padding’, ‘same’, ‘Name’, ‘CNN1’)
batchNormalizationLayer(‘Name’, ‘BN1’)
reluLayer(‘Name’, ‘ReLU1’)
maxPooling2dLayer(poolSize, ‘Stride’, [1 2], ‘Name’, ‘MaxPool1’)

convolution2dLayer(filterSize, 24*netWidth, ‘Padding’, ‘same’, ‘Name’, ‘CNN2’)
batchNormalizationLayer(‘Name’, ‘BN2’)
reluLayer(‘Name’, ‘ReLU2’)
maxPooling2dLayer(poolSize, ‘Stride’, [1 2], ‘Name’, ‘MaxPool2’)

convolution2dLayer(filterSize, 32*netWidth, ‘Padding’, ‘same’, ‘Name’, ‘CNN3’)
batchNormalizationLayer(‘Name’, ‘BN3’)
reluLayer(‘Name’, ‘ReLU3’)
maxPooling2dLayer(poolSize, ‘Stride’, [1 2], ‘Name’, ‘MaxPool3’)

convolution2dLayer(filterSize, 48*netWidth, ‘Padding’, ‘same’, ‘Name’, ‘CNN4’)
batchNormalizationLayer(‘Name’, ‘BN4’)
reluLayer(‘Name’, ‘ReLU4’)
maxPooling2dLayer(poolSize, ‘Stride’, [1 2], ‘Name’, ‘MaxPool4’)

convolution2dLayer(filterSize, 64*netWidth, ‘Padding’, ‘same’, ‘Name’, ‘CNN5’)
batchNormalizationLayer(‘Name’, ‘BN5’)
reluLayer(‘Name’, ‘ReLU5’)
maxPooling2dLayer(poolSize, ‘Stride’, [1 2], ‘Name’, ‘MaxPool5’)

convolution2dLayer(filterSize, 96*netWidth, ‘Padding’, ‘same’, ‘Name’, ‘CNN6’)
batchNormalizationLayer(‘Name’, ‘BN6’)
reluLayer(‘Name’, ‘ReLU6’)

averagePooling2dLayer([1 ceil(spf/32)], ‘Name’, ‘AP1’)

fullyConnectedLayer(numModTypes, ‘Name’, ‘FC1’)
softmaxLayer(‘Name’, ‘SoftMax’)

classificationLayer(‘Name’, ‘Output’) ]

% 接下来配置 TrainingOptionsSGDM 以使用小批量大小为 256 的 SGDM 求解器。将最大轮数设置为 12，因为更多轮数不会提供进一步的训练优势。通过将执行环境设置为 ‘gpu’，在 GPU 上训练网络。将初始学习率设置为 2x10?2。每 9 轮后将学习率降低十分之一。将 ‘Plots’ 设置为“training-progress’ 以对训练进度绘图。在 NVIDIA Titan Xp GPU 上，网络需要大约 25 分钟来完成训练。
maxEpochs = 12;
miniBatchSize = 256;
validationFrequency = floor(numel(rxTrainingLabel)/miniBatchSize);
options = trainingOptions(‘sgdm’, …
‘InitialLearnRate’,2e-2, …
‘MaxEpochs’,maxEpochs, …
‘MiniBatchSize’,miniBatchSize, …
‘Shuffle’,‘every-epoch’, …
‘Plots’,‘training-progress’, …
‘Verbose’,false, …
‘ValidationData’,{rxValidation,rxValidationLabel}, …
‘ValidationFrequency’,validationFrequency, …
‘LearnRateSchedule’, ‘piecewise’, …
‘LearnRateDropPeriod’, 9, …
‘LearnRateDropFactor’, 0.1, …
‘ExecutionEnvironment’, ‘gpu’);

% 或者训练网络，或者使用已经过训练的网络。默认情况下，此示例使用经过训练的网络。
if trainNow == true
fprintf(‘%s - Training the network\n’, datestr(toc/86400,‘HH:MM:SS’))
trainedNet = trainNetwork(rxTraining,rxTrainingLabel,modClassNet,options);
else
load trainedModulationClassificationNetwork
end
% 如训练进度图所示，网络在大约 12 轮后收敛于几乎 90% 的准确度。

% 通过获得测试帧的分类准确度来评估经过训练的网络。结果表明，该网络对这组波形实现的准确度达到 95% 左右。
fprintf(‘%s - Classifying test frames\n’, datestr(toc/86400,‘HH:MM:SS’))
rxTestPred = classify(trainedNet,rxTest);
testAccuracy = mean(rxTestPred == rxTestLabel);
disp("Test accuracy: " + testAccuracy*100 + “%”)

% 绘制测试帧的混淆矩阵。如矩阵所示，网络混淆了 16-QAM 和 64-QAM 帧。此问题是预料之中的，因为每个帧只携带 128 个符号，而 16-QAM 是 64-QAM 的子集。该网络还混淆了 QPSK 和 8-PSK 帧，因为在信道衰落和频率偏移引发相位旋转后，这些调制类型的星座图看起来相似。
figure
cm = confusionchart(rxTestLabel, rxTestPred);
cm.Title = ‘Confusion Matrix for Test Data’;
cm.RowSummary = ‘row-normalized’;
cm.Parent.Position = [cm.Parent.Position(1:2) 740 424];

%% 3 进一步探查
% 要提高准确度，可以优化网络参数，例如过滤器数量、过滤器大小；或者优化网络结构，例如添加更多层，使用不同激活层等。
% Communication Toolbox 提供了更多调制类型和信道减损。有关详细信息，请参阅Modulation (Communications Toolbox) 和部分。您还可以使用 LTE Toolbox、WLAN Toolbox 和 5G Toolbox 添加标准特定信号。您还可以使用 Phased Array System Toolbox 添加雷达信号。
% 附录：“调制器”部分提供用于生成调制信号的 MATLAB 函数。您还可以探查以下函数和 System object 以获得详细信息：
% helperModClassTestChannel.m
% helperModClassFrameGenerator.m
% helperModClassFrameStore.m
% 4 参考文献
% O’Shea, T. J., J. Corgan, and T. C. Clancy. “Convolutional Radio Modulation Recognition Networks.” Preprint, submitted June 10, 2016. https://arxiv.org/abs/1602.04105
% O’Shea, T. J., T. Roy, and T. C. Clancy. “Over-the-Air Deep Learning Based Radio Signal Classification.” IEEE Journal of Selected Topics in Signal Processing. Vol. 12, Number 1, 2018, pp. 168–179.
% Liu, X., D. Yang, and A. E. Gamal. “Deep Neural Network Architectures for Modulation Classification.” Preprint, submitted January 5, 2018. https://arxiv.org/abs/1712.00443v3

%% 5 附录：辅助函数
function modulator = getModulator(modType, sps, fs)
%getModulator Modulation function selector
% MOD = getModulator(TYPE,SPS,FS) returns the modulator function handle
% MOD based on TYPE. SPS is the number of samples per symbol and FS is
% the sample rate.

switch modType
case “BPSK”
modulator = @(x)bpskModulator(x,sps);
case “QPSK”
modulator = @(x)qpskModulator(x,sps);
case “8PSK”
modulator = @(x)psk8Modulator(x,sps);
case “16QAM”
modulator = @(x)qam16Modulator(x,sps);
case “64QAM”
modulator = @(x)qam64Modulator(x,sps);
case “GFSK”
modulator = @(x)gfskModulator(x,sps);
case “CPFSK”
modulator = @(x)cpfskModulator(x,sps);
case “PAM4”
modulator = @(x)pam4Modulator(x,sps);
case “B-FM”
modulator = @(x)bfmModulator(x, fs);
case “DSB-AM”
modulator = @(x)dsbamModulator(x, fs);
case “SSB-AM”
modulator = @(x)ssbamModulator(x, fs);
end
end

⛄三、运行结果

⛄四、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1]桂祥胜,洪居亭,代华建,孙田亮.一种基于卷积神经网络的信号调制方式识别方法[J].现代计算机(专业版). 2019,(10)

3 备注
简介此部分摘自互联网，仅供参考，若侵权，联系删除