搜索關(guān)于信噪比與向已知信號添加噪聲的整理及個人心得
2017-03-19 by:CAE仿真在線 來源:互聯(lián)網(wǎng)
首先聲明,我不是專門搞這玩意的,只不過最近要在信號中添加不同的噪聲,所以收集了一些資料,也積累了一些心得。歡迎大家指正。
說起“向已知信號添加噪聲”,有一個帖子不得不提,那是由會員happy介紹的兩個常用函數(shù),我這里引用一下:
MATLAB中產(chǎn)生高斯白噪聲非常方便,可以直接應用兩個函數(shù),一個是WGN,另一個是AWGN。WGN用于產(chǎn)生高斯白噪聲,AWGN則用于在某一信號中加入高斯白噪聲。
1. WGN:產(chǎn)生高斯白噪聲
y = wgn(m,n,p) 產(chǎn)生一個m行n列的高斯白噪聲的矩陣,p以dBW為單位指定輸出噪聲的強度。
y = wgn(m,n,p,imp) 以歐姆(Ohm)為單位指定負載阻抗。
y = wgn(m,n,p,imp,state) 重置RANDN的狀態(tài)。
在數(shù)值變量后還可附加一些標志性參數(shù):
y = wgn(…,POWERTYPE) 指定p的單位。POWERTYPE可以是'dBW', 'dBm'或'linear'。線性強度(linear power)以瓦特(Watt)為單位。
y = wgn(…,OUTPUTTYPE) 指定輸出類型。OUTPUTTYPE可以是'real'或'complex'。
2. AWGN:在某一信號中加入高斯白噪聲
y = awgn(x,SNR) 在信號x中加入高斯白噪聲。信噪比SNR以dB為單位。x的強度假定為0dBW。如果x是復數(shù),就加入復噪聲。
y = awgn(x,SNR,SIGPOWER) 如果SIGPOWER是數(shù)值,則其代表以dBW為單位的信號強度;如果SIGPOWER為'measured',則函數(shù)將在加入噪聲之前測定信號強度。
y = awgn(x,SNR,SIGPOWER,STATE) 重置RANDN的狀態(tài)。
y = awgn(…,POWERTYPE) 指定SNR和SIGPOWER的單位。POWERTYPE可以是'dB'或'linear'。如果POWERTYPE是'dB',那么SNR以dB為單位,而SIGPOWER以dBW為單位。如果POWERTYPE是'linear',那么SNR作為比值來度量,而SIGPOWER以瓦特為單位。
注釋:
1. 分貝(decibel, dB):分貝(dB)是表示相對功率或幅度電平的標準單位,換句話說,就是我們用來表示兩個能量之間的差別的一種表示單位,它不是一個絕對單位。例如,電子系統(tǒng)中將電壓、電流、功率等物理量的強弱通稱為電平,電平的單位通常就以分貝表示,即事先取一個電壓或電流作為參考值(0dB),用待表示的量與參考值之比取對數(shù),再乘以20作為電平的分貝數(shù)(功率的電平值改乘10)。
2. 分貝瓦(dBW, dB Watt):指以1W的輸出功率為基準時,用分貝來測量的功率放大器的功率值。
3. dBm (dB-milliWatt):即與1milliWatt(毫瓦)作比較得出的數(shù)字。
0 dBm = 1 mW
10 dBm = 10 mW
20 dBm = 100 mW
也可直接用randn函數(shù)產(chǎn)生高斯分布序列,例如:
y=randn(1,2500);
y=y/std(y);
y=y-mean(y);
a=0.0128;
b=sqrt(0.9596);
y=a+b*y;
就得到了 N ( 0.0128, 0.9596 ) 的高斯分布序列
產(chǎn)生指定方差和均值的隨機數(shù)
設某個隨機變量x均值為mu,方差為var^2,若要產(chǎn)生同樣分布的隨機變量y,但使新的隨機變量參數(shù)改變,均值為mu_1,方差為var_1^2,可以用如下公式進行變換:
y=var_1/var*(x-mu)+mu_1,其中x為隨機變量,其余為常數(shù)(原分布參數(shù))。
具體到正態(tài)分布,若要產(chǎn)生均值為u,方差為o^2的M*N的隨機數(shù)矩陣,可以用
y=o*randn(M,N)+u得到。
對于均勻分布,若要產(chǎn)生[a,b]區(qū)間的均勻分布的M*N的隨機數(shù)矩陣,則可以用
y=rand(M,N)*(b-a)+a得到。
上述資料基本上完整地描述了原始問題,不過有幾點內(nèi)容附帶說明一下:
1. 首先更正一個錯誤,我認為在“生成N ( 0.0128, 0.9596 ) 的高斯分布序列”的程序中,應該改為以下的代碼:
y=randn(1,2500);
y=y-mean(y);
y=y/std(y);
a=0.0128;
b=sqrt(0.9596);
y=a+b*y;
注:經(jīng)驗證,“先除方差后減均值” 與 “先減均值后除方差” 兩者是一致的
2. 上面資料最后部分隱含了一個出自zhyuer 版友的結(jié)論:
1) rand產(chǎn)生的是[0,1]上的均勻分布的隨機序列
2) randn產(chǎn)生均值為0,方差為1的高斯隨機序列,也就是白噪聲序列;
也就是說,可以直接使用上面兩個函數(shù)對原始信號添加噪聲(例如y=x+rand(length(x),1)或者y=x+randn(length(x),1))
3. 事實上,無論是wgn還是awgn函數(shù),實質(zhì)都是由randn函數(shù)產(chǎn)生的噪聲。即,wgn函數(shù)中調(diào)用了randn函數(shù),而awgn函數(shù)中調(diào)用了wgn函數(shù)。下面就我熟悉的“向已知信號添加某個信噪比(SNR)的高斯白噪聲”來說明一下,不過如果大家閱讀過awgn的實現(xiàn)代碼就不用看下去了,呵呵。
從上述可知,這個任務可以使用awgn函數(shù)實現(xiàn),具體命令是:awgn(x,snr,’measured’,'linear’),命令的作用是對原信號f(x)添加信噪比(比值)為SNR的噪聲,在添加之前先估計信號f的強度。這里涉及三個問題:在awgn這個函數(shù)中,SNR是如何計算的?什么是信號的強度?awgn函數(shù)具體是如何添加噪聲的?
事實上,前兩個問題是相關(guān)的,因為根據(jù)定義,SNR就是信號的強度除以噪聲的強度,所以,首先來講講信號的強度。其實信號的強度指的就是信號的能量,在連續(xù)的情形就是對f(x)平方后求積分,而在離散的情形自然是求和代替積分了。在matlab中也是這樣實現(xiàn)的,只不過多了一個規(guī)范化步驟罷了:
sigPower = sum(abs(sig(:)).^2)/length(sig(:))
這就是信號的強度。至此,SNR的具體實現(xiàn)也不用多說了(注:由于采用的是比值而非db,所以與下面“計算信噪比”所使用的方式不同,即沒有求對數(shù)步驟)。
最后說說awgn函數(shù)具體是如何添加噪聲的。
事實上也很簡單,在求出f的強度后,結(jié)合指定的信噪比,就可以求出需要添加的噪聲的強度noisePower=sigPower/SNR。由于使用的是高斯白噪聲即randn函數(shù),而randn的結(jié)果是一個強度為1的隨機序列(自己試試sum(randn(1000,1).^2)/1000就知道了,注意信號的長度不能太小)。于是,所要添加的噪聲信號顯然就是:sqrt(noisePower)*randn(n,1),其中n為信號長度。
4. 上面所說的都是具有分布特性(相關(guān)的)隨機序列,如果需要添加不相關(guān)的隨機序列,則可以使用jimin版友的方法:
for i=1:100
x(i)=randn(1);
end
即先產(chǎn)生噪聲信號,后再與原信號疊加。
最后是另外的一些常見問題,整理如下:
1. Matlab中如何產(chǎn)生值為0,1的隨機序列?【轉(zhuǎn)bainhome版友】:round(rand(5))
2. Matlab中如何計算信噪比?下面的代碼轉(zhuǎn)自會員Happy:
function snr=SNR(I,In)
% 計算信號噪聲比函數(shù)
% by Qulei
% I :original signal
% In:noisy signal
% ie. Original signal + noise signal
% snr=10*log10(sigma2(I2)/sigma2(I2-I1))
[row,col,nchannel]=size(I);
snr=0;
if nchannel==1
%gray image
Ps=sum(sum((I-mean(mean(I))).^2));
%signal power
Pn=sum(sum((I-In).^2));
%noise power
snr=10*log10(Ps/Pn);
elseif nchannel==3
%color image
for i=1:3
Ps=sum(sum((I(:,:,i)-mean(mean(I(:,:,i)))).^2));
%signal power
Pn=sum(sum((I(:,:,i)-In(:,:,i)).^2));
%noise power
snr=snr+10*log10(Ps/Pn);
end
snr=snr/3;
end
3. 隨機產(chǎn)生1-n的索引排列:randperm函數(shù)
4. 隨機產(chǎn)生1-60的正整數(shù)一個(三種方法):
A=randperm(60);
b=A(1)
A=round(100*rand(1,10))
b=A(1)
b = unidrnd(60,1,1)
聲明:本文由聲振論壇會員eight原創(chuàng),版權(quán)歸原作者所有,轉(zhuǎn)載請聯(lián)系原作者授權(quán),并注明:來自@聲振之家。
相關(guān)標簽搜索：關(guān)于信噪比與向已知信號添加噪聲的整理及個人心得 Actran培訓 Actran培訓課程 Actran分析 Actran在線視頻教程 Actran技術(shù)學習教程 Actran軟件教程 Actran資料下載 Actran代做 Actran基礎(chǔ)知識 Fluent、CFX流體分析 HFSS電磁分析 Ansys培訓
