Классификация модуляций с глубоким обучением

Этот пример использует:

В этом примере показано, как использовать сверточную нейронную сеть (CNN) в классификации модуляций. Вы генерируете синтетический продукт, поврежденные каналом формы волны. Используя сгенерированные формы волны как обучающие данные, вы обучаете CNN классификации модуляций. Вы затем тестируете CNN с оборудованием программно определяемого радио (SDR) и беспроводными сигналами.

Предскажите тип модуляции Используя CNN

Обученный CNN в этом примере распознает эти цифровые восемь и три аналоговых типа модуляции:

Бинарное манипулирование сдвига фазы (BPSK)
Манипулирование сдвига фазы Quadrature (QPSK)
8-ary сдвиг фазы, включающий (8-PSK)
16-ary квадратурная амплитудная (16-QAM) модуляция
64-ary квадратурная амплитудная (64-QAM) модуляция
4-ary импульсная амплитудная модуляция (PAM4)
Гауссово манипулирование сдвига частоты (GFSK)
Непрерывное манипулирование сдвига частоты фазы (CPFSK)
Широковещательно передайте FM (B-FM)
Двойная амплитудная модуляция боковой полосы (DSB-AM)
Однополосная амплитудная модуляция (SSB-AM)

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

Во-первых, загрузите обучивший сеть. Для получения дополнительной информации на сетевом обучении, смотрите раздел Training a CNN.

load trainedModulationClassificationNetwork
trainedNet

trainedNet = 
  SeriesNetwork with properties:

         Layers: [28x1 nnet.cnn.layer.Layer]
     InputNames: {'Input Layer'}
    OutputNames: {'Output'}

Обученный CNN берет 1 024 поврежденных каналом выборки и предсказывает тип модуляции каждой системы координат. Сгенерируйте несколько систем координат BPSK, которым повреждают с Rician многопутевое исчезновение, центральная частота и дрейф времени выборки и AWGN. Используйте randi функция, чтобы сгенерировать случайные биты, pskmod функционируйте к BPSK - модулируют биты, rcosdesign функционируйте, чтобы спроектировать повышенный формирующий фильтр импульса косинуса квадратного корня и filter функционируйте к форме импульса символы. Затем используйте CNN, чтобы предсказать тип модуляции систем координат.

% Set the random number generator to a known state to be able to regenerate
% the same frames every time the simulation is run
rng(123456)
% Random bits
d = randi([0 1],1024,1);
% BPSK modulation
syms = pskmod(d,2);
% Square-root raised cosine filter
filterCoeffs = rcosdesign(0.35,4,8);
tx = filter(filterCoeffs,1,upsample(syms,8));
% Channel
channel = helperModClassTestChannel(...
  'SampleRate',200e3, ...
  'SNR',30, ...
  'PathDelays',[0 1.8 3.4] / 200e3, ...
  'AveragePathGains',[0 -2 -10], ...
  'KFactor',4, ...
  'MaximumDopplerShift',4, ...
  'MaximumClockOffset',5, ...
  'CenterFrequency',902e6);
rx = channel(tx);
% Plot transmitted and received signals
scope = dsp.TimeScope(2,200e3,'YLimits',[-1 1],'ShowGrid',true,...
  'LayoutDimensions',[2 1],'TimeSpan',45e-3);
scope(tx,rx)

% Frame generation for classification
unknownFrames = getNNFrames(rx,'Unknown');
% Classification
[prediction1,score1] = classify(trainedNet,unknownFrames);

Возвратите прогнозы классификатора, которые походят на трудные решения. Сеть правильно идентифицирует системы координат как системы координат BPSK. Для получения дополнительной информации на генерации модулируемых сигналов, см. Приложение: Модуляторы.

prediction1

prediction1 = 7x1 categorical array
     BPSK 
     BPSK 
     BPSK 
     BPSK 
     BPSK 
     BPSK 
     BPSK

Классификатор также возвращает вектор музыки к каждой системе координат. Счет соответствует вероятности, что каждая система координат имеет предсказанный тип модуляции. Постройте баллы.

plotScores(score1,modulationTypes)

Затем используйте CNN, чтобы классифицировать системы координат PAM4.

% Random bits
d = randi([0 3], 1024, 1);
% PAM4 modulation
syms = pammod(d,4);
% Square-root raised cosine filter
filterCoeffs = rcosdesign(0.35, 4, 8);
tx = filter(filterCoeffs, 1, upsample(syms,8));
% Channel
rx = channel(tx);
% Plot transmitted and received signals
scope = dsp.TimeScope(2,200e3,'YLimits',[-2 2],'ShowGrid',true,...
  'LayoutDimensions',[2 1],'TimeSpan',45e-3);
scope(tx,rx)

% Frame generation for classification
unknownFrames = getNNFrames(rx,'Unknown');
% Classification
[estimate2,score2] = classify(trainedNet,unknownFrames);
estimate2

estimate2 = 7x1 categorical array
     PAM4 
     PAM4 
     PAM4 
     PAM4 
     PAM4 
     PAM4 
     PAM4

plotScores(score2,modulationTypes)

Прежде чем мы сможем использовать CNN в классификации модуляций или любую другую задачу, мы сначала должны обучить сеть с известным (или помеченный) данные. Первая часть этого примера показывает, как использовать функции Communications Toolbox, такие как модуляторы, фильтры, и нарушения канала, чтобы сгенерировать синтетические обучающие данные. Вторая часть фокусируется на определении, обучении и тестировании CNN для задачи классификации модуляций. Третья часть проверяет производительность сети с беспроводными сигналами с помощью платформ программно определяемого радио (SDR).

Генерация сигналов для обучения

Сгенерируйте 10 000 систем координат для каждого типа модуляции, где 80% используются в обучении, 10% используется в валидации, и 10% используется в тестировании. Мы используем системы координат обучения и валидации во время сетевой учебной фазы. Итоговая точность классификации получена с помощью тестовых систем координат. Каждая система координат является 1 024 выборками долго и имеет частоту дискретизации 200 кГц. Для цифровых типов модуляции восемь выборок представляют символ. Сеть принимает каждое решение на основе одной систем координат, а не на нескольких последовательных системах координат (как в видео). Примите центральную частоту 902 МГц и 100 МГц для цифровых и аналоговых типов модуляции, соответственно.

Чтобы запустить этот пример быстро, используйте обучивший сеть и сгенерируйте небольшое количество учебных систем координат. Чтобы обучить сеть на вашем компьютере, выберите опцию "Train network now" (т.е. установите trainNow на истину).

trainNow = false;
if trainNow == верный
  numFramesPerModType = 10000;
else
  numFramesPerModType = 500;
end
percentTrainingSamples = 80;
percentValidationSamples = 10;
percentTestSamples = 10;

SPS = 8;                % Samples per symbol
солнцезащитный фактор = 1024;             % Samples per frame
symbolsPerFrame = солнцезащитный фактор / SPS;
фс = 200e3;             % Sample rate
ФК = [902e6 100e6];     % Center frequencies

Создайте нарушения канала

Передайте каждый кадр через канал с

AWGN
Rician многопутевое исчезновение
Смещение часов, приводящее к центральному смещению частоты и дрейфу времени выборки

Поскольку сеть в этом примере принимает решения на основе одной систем координат, каждая система координат должна пройти через независимый канал.

AWGN

Канал добавляет AWGN с ОСШ 30 дБ. Поскольку системы координат нормированы, шумовое стандартное отклонение, может быть вычислен как

SNR = 30;
std = sqrt(10.^(-SNR/10))

std = 0.0316

Реализуйте канал с помощью comm.AWGNChannel,

awgnChannel = comm.AWGNChannel(...
  'NoiseMethod', 'Signal to noise ratio (SNR)', ...
  'SignalPower', 1, ...
  'SNR', SNR)

awgnChannel = 
  comm.AWGNChannel with properties:

     NoiseMethod: 'Signal to noise ratio (SNR)'
             SNR: 30
     SignalPower: 1
    RandomStream: 'Global stream'

Многопутевой Rician

Канал передает сигналы через Rician многопутевой исчезающий канал с помощью comm.RicianChannel Системный объект. Примите профиль задержки [0 1.8 3.4] выборки с соответствующими средними усилениями пути [0 - 2 - 10] дБ. K-фактор равняется 4, и максимальный эффект Доплера составляет 4 Гц, который эквивалентен скорости обхода на уровне 902 МГц. Реализуйте канал со следующими настройками.

multipathChannel = comm.RicianChannel(...
  'SampleRate', fs, ...
  'PathDelays', [0 1.8 3.4]/fs, ...
  'AveragePathGains', [0 -2 -10], ...
  'KFactor', 4, ...
  'MaximumDopplerShift', 4)

multipathChannel = 
  comm.RicianChannel with properties:

                SampleRate: 200000
                PathDelays: [0 9.0000e-06 1.7000e-05]
          AveragePathGains: [0 -2 -10]
        NormalizePathGains: true
                   KFactor: 4
    DirectPathDopplerShift: 0
    DirectPathInitialPhase: 0
       MaximumDopplerShift: 4
           DopplerSpectrum: [1x1 struct]

  Show all properties

Смещение часов

Смещение часов происходит из-за погрешностей внутренних источников часов передатчиков и получателей. Синхронизируйте причины смещения центральная частота, которая используется к downconvert сигнал к основной полосе и уровень выборки цифро-аналогового преобразователя, чтобы отличаться от идеальных значений. Средство моделирования канала использует фактор смещения часов $C$ , выраженный как $C = 1 + \frac{Δ_{часы}}{10^{6}}$ , где $Δ_{часы}$ смещение часов. Для каждой системы координат канал генерирует случайное $Δ_{часы}$ значение от равномерно распределенного множества значений в области значений [ ${- \max Δ}_{часы}$ ${\max Δ}_{часы}$ ], где ${\max Δ}_{часы}$ максимальное смещение часов. Смещение часов измеряется в частях на миллион (ppm). В данном примере примите максимальное смещение часов 5 страниц в минуту.

maxDeltaOff = 5;
deltaOff = (rand()*2*maxDeltaOff) - maxDeltaOff;
C = 1 + (deltaOff/1e6);

Смещение частоты

Подвергните каждую систему координат смещению частоты на основе фактора смещения часов $C$ и центральная частота. Реализованный канал с помощью comm.PhaseFrequencyOffset.

offset = -(C-1)*fc(1);
frequencyShifter = comm.PhaseFrequencyOffset(...
  'SampleRate', fs, ...
  'FrequencyOffset', offset)

frequencyShifter = 
  comm.PhaseFrequencyOffset with properties:

              PhaseOffset: 0
    FrequencyOffsetSource: 'Property'
          FrequencyOffset: -2.4386e+03
               SampleRate: 200000

Выборка смещения уровня

Подвергните каждую систему координат смещению уровня выборки на основе фактора смещения часов $C$ . Реализуйте канал с помощью interp1 функция, чтобы передискретизировать систему координат на новом уровне $C \times f_{s}$ .

Объединенный канал

Используйте объект 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)

channel = 
  helperModClassTestChannel with properties:

                    SNR: 30
        CenterFrequency: 902000000
             SampleRate: 200000
             PathDelays: [0 9.0000e-06 1.7000e-05]
       AveragePathGains: [0 -2 -10]
                KFactor: 4
    MaximumDopplerShift: 4
     MaximumClockOffset: 5

Можно просмотреть основную информацию о канале с помощью функции объекта info.

chInfo = info(channel)

chInfo = struct with fields:
               ChannelDelay: 6
     MaximumFrequencyOffset: 4510
    MaximumSampleRateOffset: 1

Генерация сигналов

Создайте цикл, который генерирует поврежденные каналом системы координат для каждого типа модуляции и хранит системы координат их соответствующими метками в frameStore. Удалите случайное число выборок с начала каждой системы координат удалить переходные процессы и убедиться, что системы координат имеют случайную начальную точку относительно контуров символа.

% Set the random number generator to a known state to be able to regenerate
% the same frames every time the simulation is run
rng(1235)
tic

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

00:00:00 - Generating BPSK frames
00:00:03 - Generating QPSK frames
00:00:05 - Generating 8PSK frames
00:00:08 - Generating 16QAM frames
00:00:11 - Generating 64QAM frames
00:00:14 - Generating PAM4 frames
00:00:17 - Generating GFSK frames
00:00:19 - Generating CPFSK frames
00:00:22 - Generating B-FM frames
00:00:37 - Generating DSB-AM frames
00:00:39 - Generating SSB-AM frames

Затем разделите системы координат на обучение, валидацию и тестовые данные. По умолчанию, frameStore места основополосные выборки I/Q в строках в выходных системах координат. Выходные системы координат имеют размер [2xspfx1xN], где первая строка является синфазными выборками, и вторая строка является квадратурными выборками.

[mcfsTraining,mcfsValidation,mcfsTest] = splitData(frameStore,...
  [percentTrainingSamples,percentValidationSamples,percentTestSamples]);
[rxTraining,rxTrainingLabel] = get(mcfsTraining);
[rxValidation,rxValidationLabel] = get(mcfsValidation);
[rxTest,rxTestLabel] = get(mcfsTest);

% 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")

Обучите CNN

Этот пример использует CNN, который состоит из шести слоев свертки и одного полносвязного слоя. Каждый слой свертки кроме последнего сопровождается пакетным слоем нормализации, исправил линейный модуль (ReLU) слой активации и макс. объединение слоя. В последнем слое свертки макс. слой объединения заменяется средним слоем объединения. Выходной слой имеет softmax активацию. Для руководства проектирования сети смотрите Советы Глубокого обучения и Приемы (Deep Learning Toolbox).

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') ]

modClassNet = 
  28x1 Layer array with layers:

     1   'Input Layer'   Image Input             2x1024x1 images
     2   'CNN1'          Convolution             16 1x8 convolutions with stride [1  1] and padding 'same'
     3   'BN1'           Batch Normalization     Batch normalization
     4   'ReLU1'         ReLU                    ReLU
     5   'MaxPool1'      Max Pooling             1x2 max pooling with stride [1  2] and padding [0  0  0  0]
     6   'CNN2'          Convolution             24 1x8 convolutions with stride [1  1] and padding 'same'
     7   'BN2'           Batch Normalization     Batch normalization
     8   'ReLU2'         ReLU                    ReLU
     9   'MaxPool2'      Max Pooling             1x2 max pooling with stride [1  2] and padding [0  0  0  0]
    10   'CNN3'          Convolution             32 1x8 convolutions with stride [1  1] and padding 'same'
    11   'BN3'           Batch Normalization     Batch normalization
    12   'ReLU3'         ReLU                    ReLU
    13   'MaxPool3'      Max Pooling             1x2 max pooling with stride [1  2] and padding [0  0  0  0]
    14   'CNN4'          Convolution             48 1x8 convolutions with stride [1  1] and padding 'same'
    15   'BN4'           Batch Normalization     Batch normalization
    16   'ReLU4'         ReLU                    ReLU
    17   'MaxPool4'      Max Pooling             1x2 max pooling with stride [1  2] and padding [0  0  0  0]
    18   'CNN5'          Convolution             64 1x8 convolutions with stride [1  1] and padding 'same'
    19   'BN5'           Batch Normalization     Batch normalization
    20   'ReLU5'         ReLU                    ReLU
    21   'MaxPool5'      Max Pooling             1x2 max pooling with stride [1  2] and padding [0  0  0  0]
    22   'CNN6'          Convolution             96 1x8 convolutions with stride [1  1] and padding 'same'
    23   'BN6'           Batch Normalization     Batch normalization
    24   'ReLU6'         ReLU                    ReLU
    25   'AP1'           Average Pooling         1x32 average pooling with stride [1  1] and padding [0  0  0  0]
    26   'FC1'           Fully Connected         11 fully connected layer
    27   'SoftMax'       Softmax                 softmax
    28   'Output'        Classification Output   crossentropyex

Используйте analyzeNetwork функционируйте, чтобы отобразить интерактивную визуализацию сетевой архитектуры, обнаружить ошибки и проблемы с сетью, и получить подробную информацию о сетевых слоях. Эта сеть имеет 98,323 learnables.

analyzeNetwork(modClassNet)

Затем сконфигурируйте TrainingOptionsSGDM использовать решатель SGDM с мини-пакетным размером 256. Определите максимальный номер эпох к 12, поскольку большее число эпох не обеспечивает дальнейшего учебного преимущества. Обучите сеть на графическом процессоре путем установки среды выполнения на 'gpu'. Установите начальный темп обучения на $2 x 10^{- 2}$ . Уменьшайте темп обучения фактором 10 каждых 9 эпох. Установите 'Plots' к 'training-progress' построить процесс обучения. На Титане NVIDIA Xp графический процессор сеть занимает приблизительно 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%-й точностью.

Оцените обучивший сеть путем получения точности классификации для тестовых систем координат. Результаты показывают, что сеть достигает приблизительно 90%-й точности для этой группы форм волны.

fprintf('%s - Classifying test frames\n', datestr(toc/86400,'HH:MM:SS'))

00:01:07 - Classifying test frames

rxTestPred = classify(trainedNet,rxTest);
testAccuracy = mean(rxTestPred == rxTestLabel);
disp("Test accuracy: " + testAccuracy*100 + "%")

Test accuracy: 90.5455%

Постройте матрицу беспорядка для тестовых систем координат. Когда матрица показывает, сеть путает 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];

I/Q как Страницы

По умолчанию, frameStore места основополосные выборки I/Q в строках в 2D массиве. Поскольку фильтры свертки имеют размер [1xsps], сверточный синфазный процесс слоев и квадратурные компоненты независимо. Только в полносвязном слое информация от синфазных компонентов и квадратурных объединенных компонентов.

Альтернатива должна представлять выборки I/Q как трехмерный массив размера [1xSPFx2], который помещает синфазные компоненты и квадратурные компоненты в 3-й размерности (страницы). Этот подход смешивает информацию во мне и Q даже в сверточных слоях и лучше использует информацию о фазе. Установите 'OutputFormat' свойство системы координат хранит к "IQAsPages", и размер входного слоя к [1xSPFx2].

% Put the data in [1xspfx2] format
mcfsTraining.OutputFormat = "IQAsPages";
[rxTraining,rxTrainingLabel] = get(mcfsTraining);
mcfsValidation.OutputFormat = "IQAsPages";
[rxValidation,rxValidationLabel] = get(mcfsValidation);
mcfsTest.OutputFormat = "IQAsPages";
[rxTest,rxTestLabel] = get(mcfsTest);

% Set the options
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');

% Set the input layer input size to [1xspfx2]
modClassNet(1) = ...
  imageInputLayer([1 spf 2], 'Normalization', 'none', 'Name', 'Input Layer');

Следующий анализ сети показывает, что каждый сверточный фильтр имеет размерность 1x8x2, который позволяет сверточному слою использовать и меня и данные Q в вычислении одного фильтра выход.

analyzeNetwork(modClassNet)

% Train or load the pretrained modified network
if trainNow == true
  fprintf('%s - Training the network\n', datestr(toc/86400,'HH:MM:SS'))
  trainedNet = trainNetwork(rxTraining,rxTrainingLabel,modClassNet,options);
else
  load trainedModulationClassificationNetwork2
end

Когда график процесса обучения показывает, сеть сходится приблизительно в 12 эпох больше чем с 95%-й точностью. При представлении компонентов I/Q, когда страницы вместо строк могут улучшить точность сети приблизительно на 5%.

Оцените обучивший сеть путем получения точности классификации для тестовых систем координат. Результаты показывают, что сеть достигает приблизительно 95%-й точности для этой группы форм волны.

fprintf('%s - Classifying test frames\n', datestr(toc/86400,'HH:MM:SS'))

00:01:14 - Classifying test frames

rxTestPred = classify(trainedNet,rxTest);
testAccuracy = mean(rxTestPred == rxTestLabel);
disp("Test accuracy: " + testAccuracy*100 + "%")

Test accuracy: 95.2727%

Постройте матрицу беспорядка для тестовых систем координат. Когда матрица показывает, представляя компоненты I/Q, когда страницы вместо строк существенно увеличивают способность сети точно дифференцировать 16-QAM и 64-QAM системы координат и системы координат 8-PSK и QPSK.

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];

Протестируйте с SDR

Проверьте производительность обучившего сеть с беспроводными сигналами с помощью функции sdrTest. Чтобы выполнить этот тест, вы, должно быть, выделили SDRs для передачи и приема. Можно использовать два радио ADALM-PLUTO, или одно радио ADALM-PLUTO для передачи и один радио USRP® для приема. Необходимо установить Пакет Поддержки Communications Toolbox для Радио ADALM-PLUTO. Если вы используете радио USRP®, необходимо также установить Пакет Поддержки Communications Toolbox для Радио USRP®. sdrTest функционируйте использует те же функции модуляции, как используется в генерации учебных сигналов, и затем передает их использующий радио ADALM-PLUTO. Вместо того, чтобы симулировать канал, получите поврежденные каналом сигналы с помощью SDR, который сконфигурирован для приема сигнала (ADALM-PLUTO или радио USRP®). Используйте обучивший сеть с тем же classify функция раньше ранее предсказывала тип модуляции. Выполнение следующего сегмента кода производит матрицу беспорядка и распечатывает тестовую точность.

radioPlatform = "ADALM-PLUTO";

switch radioPlatform
  case "ADALM-PLUTO"
    if isPlutoSDRInstalled () == верный
      радио = findPlutoRadio ();
      if длина (радио)> = 2
        sdrTest (радио);
      else
        disp'Selected radios not found. Skipping over-the-air test.')
      end
    end
  case {"USRP B2xx","USRP X3xx","USRP N2xx"}
    if (isUSRPInstalled () == верный) && (isPlutoSDRInstalled () == верный)
      txRadio = findPlutoRadio ();
      rxRadio = findsdru ();
      switch radioPlatform
        case "USRP B2xx"
          idx = содержит ({rxRadio. Платформа}, {'B200','B210'});
        case "USRP X3xx"
          idx = содержит ({rxRadio. Платформа}, {'X300','X310'});
        case "USRP N2xx"
          idx = содержит ({rxRadio. Платформа}, 'N200/N210/USRP2');
      end
      rxRadio = rxRadio (idx);
      if (длина (txRadio)> = 1) && (длина (rxRadio)> = 1)
        sdrTest (rxRadio);
      else
        disp'Selected radios not found. Skipping over-the-air test.')
      end
    end
end

При использовании двух стационарных радио ADALM-PLUTO, разделенных приблизительно на 2 фута, сеть достигает 99%-й общей точности со следующей матрицей беспорядка. Результаты будут варьироваться на основе экспериментальной настройки.

Дальнейшее исследование

Возможно оптимизировать сетевые параметры, такие как количество фильтров, размера фильтра, или оптимизировать сетевую структуру, такую как добавление большего количества слоев, использование различных слоев активации, и т.д. улучшить точность.

Коммуникационный Тулбокс обеспечивает намного больше типов модуляции и нарушений канала. Для получения дополнительной информации смотрите разделы Моделей Модуляции и Канала. Можно также добавить стандартные определенные сигналы с LTE Toolbox, WLAN Toolbox и 5G Toolbox. Можно также добавить радарные сигналы с Phased Array System Toolbox.

Приложение: раздел Modulators обеспечивает, функции MATLAB раньше генерировали модулируемые сигналы. Можно также исследовать следующие функции и Системные объекты для получения дополнительной информации:

Ссылки

О'Ши, T. J. Дж. Коргэн и Т. К. Клэнси. "Сверточные Радио-Сети Распознавания Модуляции". Предварительно распечатайте, представленный 10 июня 2016. https://arxiv.org/abs/1602.04105
О'Ши, T. J. Т. Рой и Т. К. Клэнси. "Беспроводная основанная на глубоком обучении Радио-Классификация Сигнала". Журнал IEEE Выбранных Тем в Обработке сигналов. Издание 12, Номер 1, 2018, стр 168–179.
Лю, X., Д. Янг и А. Э. Джамаль. "Архитектура Глубокой нейронной сети для Классификации Модуляций". Предварительно распечатайте, представленный 5 января 2018. https://arxiv.org/abs/1712.00443v3

Приложение: Функции помощника

function testAccuracy = sdrTest(radios)
%sdrTest Test CNN performance with over-the-air signals
%   A = sdrTest sends test frames from an ADALM-PLUTO radio, receives using
%   an ADALM-PLUTO or USRP radio, performs classification with the trained 
%   network and returns the overall classification accuracy. Transmitting 
%   radio uses transmit-repeat capability to send the same waveform repeatedly 
%   without loading the main loop.

modulationTypes = categorical(["BPSK", "QPSK", "8PSK", ...
  "16QAM", "64QAM", "PAM4", "GFSK", "CPFSK", "B-FM"]);
load trainedModulationClassificationNetwork2 trainedNet
numFramesPerModType = 100;

sps = 8;                % Samples per symbol
spf = 1024;             % Samples per frame
fs = 200e3;             % Sample rate

txRadio = sdrtx('Pluto');
txRadio.RadioID = 'usb:0';
txRadio.CenterFrequency = 902e6;
txRadio.BasebandSampleRate = fs;

if isfield(radios, 'Platform')
  radioPlatform = "USRP";
  % Configure USRP radio as the receiver. 
  rxRadio = comm.SDRuReceiver("Platform",radios(1).Platform);
  switch radios(1).Platform
    case {"B200","B210"}
      masterClockRate = 5e6;
      rxRadio.SerialNum = radios(1).SerialNum;
    case {"N200/N210/USRP2"}
      masterClockRate = 100e6;
      rxRadio.IPAddress = radios(1).IPAddress;
    case {"X300","X310"}
      masterClockRate = 120e6;
      rxRadio.IPAddress = radios(1).IPAddress;
  end
  rxRadio.MasterClockRate = masterClockRate;
  rxRadio.DecimationFactor = masterClockRate/fs;
  radioInfo = info(rxRadio);
  maximumGain = radioInfo.MaximumGain;
  minimumGain = radioInfo.MinimumGain;
else
  radioPlatform = "PlutoSDR";
  rxRadio = sdrrx('Pluto');
  rxRadio.RadioID = 'usb:1';
  rxRadio.BasebandSampleRate = fs;
  rxRadio.ShowAdvancedProperties = true;
  rxRadio.EnableQuadratureCorrection = false;
  rxRadio.GainSource = "Manual";
  maximumGain = 73;
  minimumGain = -10;
end
rxRadio.SamplesPerFrame = spf;
% Use burst mode with numFramesInBurst set to 1, so that each capture 
% (call to the receiver) will return an independent fresh frame even 
% though the radio overruns.
rxRadio.EnableBurstMode = true;
rxRadio.NumFramesInBurst = 1;
rxRadio.OutputDataType = 'single';

% Display Tx and Rx radios
txRadio
rxRadio

% Set random number generator to a known state to be able to regenerate
% the same frames every time the simulation is run
rng(1235)
tic

numModulationTypes = length(modulationTypes);
txModType = repmat(modulationTypes(1),numModulationTypes*numFramesPerModType,1);
estimatedModType = repmat(modulationTypes(1),numModulationTypes*numFramesPerModType,1);
frameCnt = 1;
for modType = 1:numModulationTypes
  fprintf('%s - Testing %s frames\n', ...
    datestr(toc/86400,'HH:MM:SS'), modulationTypes(modType))
  dataSrc = getSource(modulationTypes(modType), sps, 2*spf, fs);
  modulator = getModulator(modulationTypes(modType), sps, fs);
  if contains(char(modulationTypes(modType)), {'B-FM'})...
           && (radioPlatform == "PlutoSDR")
    % Analog modulation types use a center frequency of 100 MHz
    txRadio.CenterFrequency = 100e6;
    rxRadio.CenterFrequency = 100e6;
  else
    % Digital modulation types use a center frequency of 902 MHz
    txRadio.CenterFrequency = 902e6;
    rxRadio.CenterFrequency = 902e6;
  end
  
  disp('Starting transmitter')
  x = dataSrc();
  y = modulator(x);
  % Remove filter transients
  y = y(4*sps+1:end,1);
  maxVal = max(max(abs(real(y))), max(abs(imag(y))));
  y = y *0.8/maxVal;
  % Download waveform signal to radio and repeatedly transmit it over the air
  transmitRepeat(txRadio, complex(y));
  
  disp('Adjusting receiver gain')
  rxRadio.Gain = maximumGain;
  gainAdjusted = false;
  while ~gainAdjusted
    for p=1:20
      rx = rxRadio();
    end
    maxAmplitude = max([abs(real(rx)); abs(imag(rx))]);
    if (maxAmplitude < 0.8) || (rxRadio.Gain <= minimumGain)
      gainAdjusted = true;
    else
      rxRadio.Gain = rxRadio.Gain - 3;
    end
  end
  
  disp('Starting receiver and test')
  for p=1:numFramesPerModType
    rx = rxRadio();
    
    frameEnergy = sum(abs(rx).^2);
    rx = rx / sqrt(frameEnergy);
    reshapedRx(1,:,1,1) = real(rx);
    reshapedRx(1,:,2,1) = imag(rx);
    
    % Classify
    txModType(frameCnt) = modulationTypes(modType);
    estimatedModType(frameCnt) = classify(trainedNet, reshapedRx);
    
    frameCnt = frameCnt + 1;
    
    % Pause for some time to get an independent channel. The pause duration 
    % together with the processing time of a single loop must be greater 
    % than the channel coherence time. Assume channel coherence time is less 
    % than 0.1 seconds.
    pause(0.1)
  end
  disp('Releasing Tx radio')
  release(txRadio);
  testAccuracy = mean(txModType(1:frameCnt-1) == estimatedModType(1:frameCnt-1));
  disp("Test accuracy: " + testAccuracy*100 + "%")
end
disp('Releasing Rx radio')
release(rxRadio);
testAccuracy = mean(txModType == estimatedModType);
disp("Final test accuracy: " + testAccuracy*100 + "%")

figure
cm = confusionchart(txModType, estimatedModType);
cm.Title = 'Confusion Matrix for Test Data';
cm.RowSummary = 'row-normalized';
cm.Parent.Position = [cm.Parent.Position(1:2) 740 424];
end

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

function src = getSource(modType, sps, spf, fs)
%getSource Source selector for modulation types
%    SRC = getSource(TYPE,SPS,SPF,FS) returns the data source
%    for the modulation type TYPE, with the number of samples
%    per symbol SPS, the number of samples per frame SPF, and
%    the sampling frequency FS.

switch modType
  case {"BPSK","GFSK","CPFSK"}
    M = 2;
    src = @()randi([0 M-1],spf/sps,1);
  case {"QPSK","PAM4"}
    M = 4;
    src = @()randi([0 M-1],spf/sps,1);
  case "8PSK"
    M = 8;
    src = @()randi([0 M-1],spf/sps,1);
  case "16QAM"
    M = 16;
    src = @()randi([0 M-1],spf/sps,1);
  case "64QAM"
    M = 64;
    src = @()randi([0 M-1],spf/sps,1);
  case {"B-FM","DSB-AM","SSB-AM"}
    src = @()getAudio(spf,fs);
end
end

function x = getAudio(spf,fs)
%getAudio Audio source for analog modulation types
%    A = getAudio(SPF,FS) returns the audio source A, with the
%    number of samples per frame SPF, and the sample rate FS.

persistent audioSrc audioRC

if isempty(audioSrc)
  audioSrc = dsp.AudioFileReader('audio_mix_441.wav',...
    'SamplesPerFrame',spf,'PlayCount',inf);
  audioRC = dsp.SampleRateConverter('Bandwidth',30e3,...
    'InputSampleRate',audioSrc.SampleRate,...
    'OutputSampleRate',fs);
  [~,decimFactor] = getRateChangeFactors(audioRC);
  audioSrc.SamplesPerFrame = ceil(spf / fs * audioSrc.SampleRate / decimFactor) * decimFactor;
end

x = audioRC(audioSrc());
x = x(1:spf,1);
end

function frames = getNNFrames(rx,modType)
%getNNFrames Generate formatted frames for neural networks
%   F = getNNFrames(X,MODTYPE) formats the input X, into frames
%   that can be used with the neural network designed in this
%   example, and returns the frames in the output F.

frames = helperModClassFrameGenerator(rx,1024,1024,32,8);
frameStore = helperModClassFrameStore(10,1024,categorical({modType}));
add(frameStore,frames,modType);
frames = get(frameStore);
end

function plotScores(score,labels)
%plotScores Plot classification scores of frames
%   plotScores(SCR,LABELS) plots the classification scores SCR as a stacked
%   bar for each frame. SCR is a matrix in which each row is the score for a
%   frame.

co = [0.08 0.9 0.49;
  0.52 0.95 0.70;
  0.36 0.53 0.96;
  0.09 0.54 0.67;
  0.48 0.99 0.26;
  0.95 0.31 0.17;
  0.52 0.85 0.95;
  0.08 0.72 0.88;
  0.12 0.45 0.69;
  0.22 0.11 0.49;
  0.65 0.54 0.71];
figure; ax = axes('ColorOrder',co,'NextPlot','replacechildren');
bar(ax,[score; nan(2,11)],'stacked'); legend(categories(labels),'Location','best');
xlabel('Frame Number'); ylabel('Score'); title('Classification Scores')
end

function plotTimeDomain(rxTest,rxTestLabel,modulationTypes,fs)
%plotTimeDomain Time domain plots of frames

numRows = ceil(length(modulationTypes) / 4);
spf = size(rxTest,2);
t = 1000*(0:spf-1)/fs;
if size(rxTest,1) == 2
  IQAsRows = true;
else
  IQAsRows = false;
end
for modType=1:length(modulationTypes)
  subplot(numRows, 4, modType);
  idxOut = find(rxTestLabel == modulationTypes(modType), 1);
  if IQAsRows
    rxI = rxTest(1,:,1,idxOut);
    rxQ = rxTest(2,:,1,idxOut);
  else
    rxI = rxTest(1,:,1,idxOut);
    rxQ = rxTest(1,:,2,idxOut);
  end
  plot(t,squeeze(rxI), '-'); grid on; axis equal; axis square
  hold on
  plot(t,squeeze(rxQ), '-'); grid on; axis equal; axis square
  hold off
  title(string(modulationTypes(modType)));
  xlabel('Time (ms)'); ylabel('Amplitude')
end
end

function plotSpectrogram(rxTest,rxTestLabel,modulationTypes,fs,sps)
%plotSpectrogram Spectrogram of frames

if size(rxTest,1) == 2
  IQAsRows = true;
else
  IQAsRows = false;
end
numRows = ceil(length(modulationTypes) / 4);
for modType=1:length(modulationTypes)
  subplot(numRows, 4, modType);
  idxOut = find(rxTestLabel == modulationTypes(modType), 1);
  if IQAsRows
    rxI = rxTest(1,:,1,idxOut);
    rxQ = rxTest(2,:,1,idxOut);
  else
    rxI = rxTest(1,:,1,idxOut);
    rxQ = rxTest(1,:,2,idxOut);
  end
  rx = squeeze(rxI) + 1i*squeeze(rxQ);
  spectrogram(rx,kaiser(sps),0,1024,fs,'centered');
  title(string(modulationTypes(modType)));
end
h = gcf; delete(findall(h.Children, 'Type', 'ColorBar'))
end

function flag = isPlutoSDRInstalled
%isPlutoSDRInstalled Check if ADALM-PLUTO Radio HSP is installed

spkg = matlabshared.supportpkg.getInstalled;
flag = ~isempty(spkg) && any(contains({spkg.Name},'ADALM-PLUTO','IgnoreCase',true));
end

function flag = isUSRPInstalled
%isUSRPInstalled Check if USRP Radio HSP is installed

spkg = matlabshared.supportpkg.getInstalled;
flag = ~isempty(spkg) && any(contains({spkg.Name},'USRP','IgnoreCase',true));
end

Приложение: модуляторы

function y = bpskModulator(x,sps)
%bpskModulator BPSK modulator with pulse shaping
%   Y = bpskModulator(X,SPS) BPSK modulates the input X, and returns the
%   root-raised cosine pulse shaped signal Y. X must be a column vector
%   of values in the set [0 1]. The root-raised cosine filter has a
%   roll-off factor of 0.35 and spans four symbols. The output signal
%   Y has unit power.

persistent filterCoeffs
if isempty(filterCoeffs)
  filterCoeffs = rcosdesign(0.35, 4, sps);
end
% Modulate
syms = pskmod(x,2);
% Pulse shape
y = filter(filterCoeffs, 1, upsample(syms,sps));
end

function y = qpskModulator(x,sps)
%qpskModulator QPSK modulator with pulse shaping
%   Y = qpskModulator(X,SPS) QPSK modulates the input X, and returns the
%   root-raised cosine pulse shaped signal Y. X must be a column vector
%   of values in the set [0 3]. The root-raised cosine filter has a
%   roll-off factor of 0.35 and spans four symbols. The output signal
%   Y has unit power.

persistent filterCoeffs
if isempty(filterCoeffs)
  filterCoeffs = rcosdesign(0.35, 4, sps);
end
% Modulate
syms = pskmod(x,4,pi/4);
% Pulse shape
y = filter(filterCoeffs, 1, upsample(syms,sps));
end

function y = psk8Modulator(x,sps)
%psk8Modulator 8-PSK modulator with pulse shaping
%   Y = psk8Modulator(X,SPS) 8-PSK modulates the input X, and returns the
%   root-raised cosine pulse shaped signal Y. X must be a column vector
%   of values in the set [0 7]. The root-raised cosine filter has a
%   roll-off factor of 0.35 and spans four symbols. The output signal
%   Y has unit power.

persistent filterCoeffs
if isempty(filterCoeffs)
  filterCoeffs = rcosdesign(0.35, 4, sps);
end
% Modulate
syms = pskmod(x,8);
% Pulse shape
y = filter(filterCoeffs, 1, upsample(syms,sps));
end

function y = qam16Modulator(x,sps)
%qam16Modulator 16-QAM modulator with pulse shaping
%   Y = qam16Modulator(X,SPS) 16-QAM modulates the input X, and returns the
%   root-raised cosine pulse shaped signal Y. X must be a column vector
%   of values in the set [0 15]. The root-raised cosine filter has a
%   roll-off factor of 0.35 and spans four symbols. The output signal
%   Y has unit power.

persistent filterCoeffs
if isempty(filterCoeffs)
  filterCoeffs = rcosdesign(0.35, 4, sps);
end
% Modulate and pulse shape
syms = qammod(x,16,'UnitAveragePower',true);
% Pulse shape
y = filter(filterCoeffs, 1, upsample(syms,sps));
end

function y = qam64Modulator(x,sps)
%qam64Modulator 64-QAM modulator with pulse shaping
%   Y = qam64Modulator(X,SPS) 64-QAM modulates the input X, and returns the
%   root-raised cosine pulse shaped signal Y. X must be a column vector
%   of values in the set [0 63]. The root-raised cosine filter has a
%   roll-off factor of 0.35 and spans four symbols. The output signal
%   Y has unit power.

persistent filterCoeffs
if isempty(filterCoeffs)
  filterCoeffs = rcosdesign(0.35, 4, sps);
end
% Modulate
syms = qammod(x,64,'UnitAveragePower',true);
% Pulse shape
y = filter(filterCoeffs, 1, upsample(syms,sps));
end

function y = pam4Modulator(x,sps)
%pam4Modulator PAM4 modulator with pulse shaping
%   Y = pam4Modulator(X,SPS) PAM4 modulates the input X, and returns the
%   root-raised cosine pulse shaped signal Y. X must be a column vector
%   of values in the set [0 3]. The root-raised cosine filter has a
%   roll-off factor of 0.35 and spans four symbols. The output signal
%   Y has unit power.

persistent filterCoeffs amp
if isempty(filterCoeffs)
  filterCoeffs = rcosdesign(0.35, 4, sps);
  amp = 1 / sqrt(mean(abs(pammod(0:3, 4)).^2));
end
% Modulate
syms = amp * pammod(x,4);
% Pulse shape
y = filter(filterCoeffs, 1, upsample(syms,sps));
end

function y = gfskModulator(x,sps)
%gfskModulator GFSK modulator
%   Y = gfskModulator(X,SPS) GFSK modulates the input X and returns the
%   signal Y. X must be a column vector of values in the set [0 1]. The
%   BT product is 0.35 and the modulation index is 1. The output signal
%   Y has unit power.

persistent mod meanM
if isempty(mod)
  M = 2;
  mod = comm.CPMModulator(...
    'ModulationOrder', M, ...
    'FrequencyPulse', 'Gaussian', ...
    'BandwidthTimeProduct', 0.35, ...
    'ModulationIndex', 1, ...
    'SamplesPerSymbol', sps);
  meanM = mean(0:M-1);
end
% Modulate
y = mod(2*(x-meanM));
end

function y = cpfskModulator(x,sps)
%cpfskModulator CPFSK modulator
%   Y = cpfskModulator(X,SPS) CPFSK modulates the input X and returns
%   the signal Y. X must be a column vector of values in the set [0 1].
%   the modulation index is 0.5. The output signal Y has unit power.

persistent mod meanM
if isempty(mod)
  M = 2;
  mod = comm.CPFSKModulator(...
    'ModulationOrder', M, ...
    'ModulationIndex', 0.5, ...
    'SamplesPerSymbol', sps);
  meanM = mean(0:M-1);
end
% Modulate
y = mod(2*(x-meanM));
end

function y = bfmModulator(x,fs)
%bfmModulator Broadcast FM modulator
%   Y = bfmModulator(X,FS) broadcast FM modulates the input X and returns
%   the signal Y at the sample rate FS. X must be a column vector of
%   audio samples at the sample rate FS. The frequency deviation is 75 kHz
%   and the pre-emphasis filter time constant is 75 microseconds.

persistent mod
if isempty(mod)
  mod = comm.FMBroadcastModulator(...
    'AudioSampleRate', fs, ...
    'SampleRate', fs);
end
y = mod(x);
end

function y = dsbamModulator(x,fs)
%dsbamModulator Double sideband AM modulator
%   Y = dsbamModulator(X,FS) double sideband AM modulates the input X and
%   returns the signal Y at the sample rate FS. X must be a column vector of
%   audio samples at the sample rate FS. The IF frequency is 50 kHz.

y = ammod(x,50e3,fs);
end

function y = ssbamModulator(x,fs)
%ssbamModulator Single sideband AM modulator
%   Y = ssbamModulator(X,FS) single sideband AM modulates the input X and
%   returns the signal Y at the sample rate FS. X must be a column vector of
%   audio samples at the sample rate FS. The IF frequency is 50 kHz.

y = ssbmod(x,50e3,fs);
end

Документация