Matlab實現加性高斯白噪聲信道（AWGN）下的digital調制格式識別分類
內容大綱
加性高斯白噪聲信道（AWGN）下的digital調制格式識別分類
(1. PSK; 2. QPSK; 3.8QAM; 4. 16QAM; 5. 32QAM; 6.64QAM)
100次獨立仿真，識別正確率 vs SNR
設計
我的實現方法是 基於高階累積量的信號特征的識別算法
調制格式識別過程如下：
信號預處理
- 去除直流成分
信號在接收端由於接收機的影響，有可能產生直流成分。直流分量在后面混頻等處理中會產生影響，因此在信號處理以前必須去除直流成分。
令s ̅(t)表示信號s(t)的均值，即
則去除直流后的信號表示為
%去除直流成分
CMAOUT = CMAOUT - mean(CMAOUT);
- 信號功率歸一化
由於信道衰落影響到接收信號的功率，提取有關幅度的特征參量時會出現不一致的情況。因此需要對接收信號進行功率歸一化，以消除信號功率的影響。令σ_x^2表示已經經過去除直流分量之后的信號x(t)的平均功率，即
那么經過功率歸一化后的信號表示為
%normalization接收信號功率歸一化
CMAOUT=CMAOUT/sqrt(mean(abs(CMAOUT).^2));
特征提取
基於高階累積量的信號處理方法，對通信信號中的加性高斯噪聲有很好的抑制能力，在低信噪比下進行信號識別也能有良好的性能，應用在信號分析領域是非常有效的。
隨機過程的k階累積量為
則根據定義，隨機過程的二階和四階累積量為
如果定義
，令
，則上面累積量的表達式化簡為：
在信號的實際處理中，要從有限的接收數據中估計信號的累積量，可以采用采樣點的平均代替理論的平均。例如，給定觀察數據r_k,k=1,2,⋯,N,則可以使用下來的估計表達式。當信號和噪聲的8階矩存在並為有限值的時候，其不同定義的4階累積量的估計是漸進無偏的一致估計。
由於處理的信號是在AWGN信道下接收的，所以對於C21這個二階累積量來說，由於是信號模的平方近似計算得到的，因此噪聲的功率會也包含進去了。所以，需要處理C21這個累積量。由於信噪比已知，則可以計算信號功率根據信噪比關系求得噪聲功率。
令snr為信噪比，則
C21的累積量的更新為
C21=C21-noise_power
s = CMAOUT;
signalpow = mean(abs(s).^2);%信號功率
noisepow = signalpow/(10^(snr(snrIndex)/10));%噪聲功率
C20_hat = mean(s.^2);
C21_hat = mean(abs(s).^2);
C21_hat = C21_hat-noisepow;%計算信號二階累積量C21時，由於C21為信號模的平方
%而我們接收的s是在AWGN信道下接收的，所以求C21時還應考慮噪聲功率。
C40_hat = mean(s.^4)-3*C20_hat^2;
%C41_hat = mean((s.^3).*conj(s))-3*C20_hat*C21_hat;
C42_hat = mean(abs(s).^4)-abs(C20_hat)^2-2*C21_hat^2;
C40_normal = C40_hat/C21_hat.^2;
%C42_normal = C42_hat/C21_hat.^2;
分類識別
從論文中得知的一些四階累積量的性質：
但是，對於實際接收到的信號，由於信號功率的影響，無法從累積量的絕對值中進行信號的區分。為了能夠各階累積量相比較，必須采取將信號功率歸一化的方法，使得不同信號的累積量具有可比性。歸一化是假設信號具有單位功率，即C_21=1,其他四階累積量利用C_21進行歸一化：
查詢可得不同調制方式的歸一化四階累積量理論值
根據以上表格內容可以做出判斷
AbsC40 = abs(C40_normal);
if(j==printJ) fprintf('%g ',AbsC40); end
if(AbsC40>=1.5)%PSK
classify = 1;
elseif(AbsC40>=0.9&&AbsC40<1.1)%QPSK
classify = 2;
elseif(AbsC40>=1.1&&AbsC40<1.3)%8QAM
classify = 3;
elseif(AbsC40>=0.67&&AbsC40<0.9)%16QAM
classify = 4;
elseif(AbsC40<0.35)%32QAM
classify = 5;
elseif(AbsC40>=0.35&&AbsC40<0.63)%64QAM
classify = 6;
end
if(classify == System.BitPerSymbol)
classify_correct = classify_correct + 1;
end
最后打印出信噪比與識別正確率的關系和對應的星座圖
classify_correct_ratio(snrIndex) = classify_correct/ClassifySetNumber*100;
end
%%繪制圖形
figure(1);subplot(2, 3, j);
plot(snr, classify_correct_ratio, '-b.');
figure(2);subplot(2, 3, j);
plot(real(CMAOUT),imag(CMAOUT),'.');
運行效果
實驗中用到的參數設置
%%參數設置
snr_mini = 5; %信噪比最小值
snr_max = 20; %信噪比最大值
TxSampleRate = 32e9; %信號的碼元速率
TxLinewidth = 0; %發射信號的載波線寬
TxCarrierRate = 0; %發射信號的載波頻率
DataSymbolNumber = 10000; %數據點的個數
ClassifySetNumber = 100; %獨立仿真的次數
printJ = 5; %需要輸出觀察的調制方式，0為不輸出
printXingZuo = 1; %是否需要打印星座圖，0為不打印
% signal generation;如果想要進行100組獨立的測試，可以建立100次循環，產生100組獨立的數據
for j = 1:6 % bit per symbol: 1. PSK; 2. QPSK; 3.8QAM; 4. 16QAM; 5. 32QAM; 6.64QAM...
System.BitPerSymbol = j;
snr = snr_mini:snr_max; %SNR信噪比的設置，單位dB
classify_correct_ratio = zeros(length(snr), 1);
for snrIndex= 1:length(snr)
if(j==printJ) fprintf('\n--------------- snr = %d ------------\n',snr(snrIndex)); end
classify_correct = 0;
for i = 1:ClassifySetNumber
Tx.SampleRate = TxSampleRate; %symbol Rate，信號的碼元速率，可以自行定義
Tx.Linewidth = TxLinewidth;%發射信號的載波的線寬，一般與信號的相位噪聲有關
Tx.Carrier = DataSymbolNumber;%發射信號的載波頻率
M = 2^System.BitPerSymbol;
從上述可以看出，一共進行了100次獨立仿真。
由於沒有找到32QAM的累積量的理論值，但是通過實驗探究可以觀察在32QAM的調制方式下。
通過輸出語句我們可以觀察到如下
當信噪比在不同的取值情況下面，由於歸一化處理，32QAM的C40的模值的理論值應為0.2
在不同的信噪比下面，識別正確率如下圖：
要求的信噪比是12—15，在信噪比的比較低的時候，識別會有一定的誤差，例如信噪比只有5的情況下比較明顯，隨着信噪比的提升，噪聲的影響減少，分類識別的方法體現出來良好的穩定性。可以看到，信噪比大於10之后的識別率幾乎都達到了100%。其中PSK的識別率是最好的。64QAM在累積量的識別方法里面低信噪比時稍微差點。
對應的調制方式星座圖如下：
其中64QAM的區分沒有特別明顯，這就是低信噪比時產生的識別困難情況。但是稍微提高信噪比就能很好的識別。
完整代碼
clear all;
clc;
%%參數設置
snr_mini = 5; %信噪比最小值
snr_max = 20; %信噪比最大值
TxSampleRate = 32e9; %信號的碼元速率
TxLinewidth = 0; %發射信號的載波線寬
TxCarrierRate = 0; %發射信號的載波頻率
DataSymbolNumber = 10000; %數據點的個數
ClassifySetNumber = 100; %獨立仿真的次數
printJ = 5; %需要輸出觀察的調制方式，0為不輸出
printXingZuo = 1; %是否需要打印星座圖，0為不打印
% signal generation;如果想要進行100組獨立的測試，可以建立100次循環，產生100組獨立的數據
for j = 1:6 % bit per symbol: 1. PSK; 2. QPSK; 3.8QAM; 4. 16QAM; 5. 32QAM; 6.64QAM...
System.BitPerSymbol = j;
snr = snr_mini:snr_max; %SNR信噪比的設置，單位dB
classify_correct_ratio = zeros(length(snr), 1);
for snrIndex= 1:length(snr)
if(j==printJ) fprintf('\n--------------- snr = %d ------------\n',snr(snrIndex)); end
classify_correct = 0;
for i = 1:ClassifySetNumber
Tx.SampleRate = TxSampleRate; %symbol Rate，信號的碼元速率，可以自行定義
Tx.Linewidth = TxLinewidth;%發射信號的載波的線寬，一般與信號的相位噪聲有關
Tx.Carrier = DataSymbolNumber;%發射信號的載波頻率
M = 2^System.BitPerSymbol;
%%信號生成
Tx.DataSymbol = randi([0 M-1],1,DataSymbolNumber);%每一次隨機產生的數據量
%數據的不同調制方式產生：這里把2^3（8QAM）的形式單獨拿出來設置，是為了實現最優的星型8QAM星座圖
if M ~= 8;
h = modem.qammod('M', M, 'SymbolOrder', 'Gray');
Tx.DataConstel = modulate(h,Tx.DataSymbol);
else
tmp = Tx.DataSymbol;
tmp2 = zeros(1,length(Tx.DataSymbol));
for kk = 1:length(Tx.DataSymbol)
switch tmp(kk)
case 0
tmp2(kk) = 1 + 1i;
case 1
tmp2(kk) = -1 + 1i;
case 2
tmp2(kk) = -1 - 1i;
case 3
tmp2(kk) = 1 - 1i;
case 4
tmp2(kk) = 1+sqrt(3);
case 5
tmp2(kk) = 0 + 1i .* (1+sqrt(3));
case 6
tmp2(kk) = 0 - 1i .* (1+sqrt(3));
case 7
tmp2(kk) = -1-sqrt(3);
end
end
Tx.DataConstel = tmp2;
clear tmp tmp2;
end
Tx.Signal = Tx.DataConstel;
%數據的載波加載，考慮到相位噪聲等
N = length(Tx.Signal);
dt = 1/Tx.SampleRate;
t = dt*(0:N-1);
Phase1 = [0, cumsum(normrnd(0,sqrt(2*pi*Tx.Linewidth/(Tx.SampleRate)), 1, N-1))];
carrier1 = exp(1i*(2*pi*t*Tx.Carrier + Phase1));
Tx.Signal = Tx.Signal.*carrier1;
Rx.Signal = awgn(Tx.Signal,snr(snrIndex),'measured');%數據在AWGN信道下的接收
%%信號識別
CMAOUT = Rx.Signal;
%去除直流成分
CMAOUT = CMAOUT - mean(CMAOUT);
%normalization接收信號功率歸一化
CMAOUT=CMAOUT/sqrt(mean(abs(CMAOUT).^2));
s = CMAOUT;
signalpow = mean(abs(s).^2);%信號功率
noisepow = signalpow/(10^(snr(snrIndex)/10));%噪聲功率
C20_hat = mean(s.^2);
C21_hat = mean(abs(s).^2);
C21_hat = C21_hat-noisepow;%計算信號二階累積量C21時，由於C21為信號模的平方
%而我們接收的s是在AWGN信道下接收的，所以求C21時還應考慮噪聲功率。
C40_hat = mean(s.^4)-3*C20_hat^2;
%C41_hat = mean((s.^3).*conj(s))-3*C20_hat*C21_hat;
C42_hat = mean(abs(s).^4)-abs(C20_hat)^2-2*C21_hat^2;
C40_normal = C40_hat/C21_hat.^2;
%C42_normal = C42_hat/C21_hat.^2;
AbsC40 = abs(C40_normal);
if(j==printJ) fprintf('%g ',AbsC40); end
if(AbsC40>=1.5)%PSK
classify = 1;
elseif(AbsC40>=0.9&&AbsC40<1.1)%QPSK
classify = 2;
elseif(AbsC40>=1.1&&AbsC40<1.3)%8QAM
classify = 3;
elseif(AbsC40>=0.67&&AbsC40<0.9)%16QAM
classify = 4;
elseif(AbsC40<0.35)%32QAM
classify = 5;
elseif(AbsC40>=0.35&&AbsC40<0.63)%64QAM
classify = 6;
end
if(classify == System.BitPerSymbol)
classify_correct = classify_correct + 1;
end
%subplot(1,7,snrIndex);%繪制原始噪聲
%plot(Rx.Signal,'.');
%plot(CMAOUT,'.');
end
classify_correct_ratio(snrIndex) = classify_correct/ClassifySetNumber*100;
end
%%繪制圖形
figure(1);subplot(2, 3, j);
plot(snr, classify_correct_ratio, '-b.');
axis([snr_mini snr_max 0 110]);
ylabel('識別正確率/%');
xlabel('信噪比/dB');
if(j == 1)
title('PSK調制方式識別');
elseif(j == 2)
title('QPSK調制方式識別');
elseif(j == 3)
title('8QAM調制方式識別');
elseif(j == 4)
title('16QAM調制方式識別');
elseif(j == 5)
title('32QAM調制方式識別');
else
title('64QAM調制方式識別');
end
if(printXingZuo==1)
figure(2);subplot(2, 3, j);
plot(real(CMAOUT),imag(CMAOUT),'.');
if(j == 1)
title('PSK調制方式星座圖');
elseif(j == 2)
title('QPSK調制方式星座圖');
elseif(j == 3)
title('8QAM調制方式星座圖');
elseif(j == 4)
title('16QAM調制方式星座圖');
elseif(j == 5)
title('32QAM調制方式星座圖');
else
title('64QAM調制方式星座圖');
end
end
end