说起“向已知信号添加噪声”,有一个帖子不得不提,那是由happy教授介绍的两个常用函数,我这里引用一下:

%=============================happy=================================%


MATLAB中产生高斯白噪声非常方便,可以直接应用两个函数,一个是WGN,另一个是AWGN。WGN用于产生高斯白噪声,AWGN则用于在某一信号中加入高斯白噪声。


1. WGN:产生高斯白噪声 
 
 y = wgn(m,n,p) 产生一个m行n列的高斯白噪声的矩阵,p以dBW为单位指定输出噪声的强度。 
 
 y = wgn(m,n,p,imp) 以欧姆(Ohm)为单位指定负载阻抗。 
 
 y = wgn(m,n,p,imp,state) 重置RANDN的状态。



在数值变量后还可附加一些标志性参数:


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是复数,就加入复噪声。


y = awgn(x,SNR,SIGPOWER) 如果SIGPOWER是数值,则其代表以dBW为单位的信号强度;如果SIGPOWER为'measured',则函数将在加入噪声之前测定信号强度。


y = awgn(x,SNR,SIGPOWER,STATE) 重置RANDN的状态。


y = awgn(…,POWERTYPE) 指定SNR和SIGPOWER的单位。POWERTYPE可以是'dB'或'linear'。如果POWERTYPE是'dB',那么SNR以dB为单位,而SIGPOWER以dBW为单位。如果POWERTYPE是'linear',那么SNR作为比值来度量,而SIGPOWER以瓦特为单位。



注释


1. 分贝(decibel, dB):分贝(dB)是表示相对功率或幅度电平的标准单位,换句话说,就是我们用来表示两个能量之间的差别的一种表示单位,它不是一个绝对单位。例如,电子系统中将电压、电流、功率等物理量的强弱通称为电平,电平的单位通常就以分贝表示,即事先取一个电压或电流作为参考值(0dB),用待表示的量与参考值之比取对数,再乘以20作为电平的分贝数(功率的电平值改乘10)。


2. 分贝瓦(dBW, dB Watt):指以1W的输出功率为基准时,用分贝来测量的功率放大器的功率值。


3. dBm (dB-milliWatt):即与1milliWatt(毫瓦)作比较得出的数字。


0 dBm = 1 mW 
 
 10 dBm = 10 mW 
 
 20 dBm = 100 mW


也可直接用randn函数产生高斯分布序列,例如:



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 ) 的高斯分布序列



产生指定方差和均值的随机数


设某个随机变量x均值为mu,方差为var^2,若要产生同样分布的随机变量y,但使新的随


机变量参数改变,均值为mu_1,方差为var_1^2,可以用如下公式进行变换:


y=var_1/var*(x-mu)+mu_1,其中x为随机变量,其余为常数(原分布参数)。


具体到正态分布,若要产生均值为u,方差为o^2的M*N的随机数矩阵,可以用


y=o*randn(M,N)+u得到。


对于均匀分布,若要产生[a,b]区间的均匀分布的M*N的随机数矩阵,则可以用


y=rand(M,N)*(b-a)+a得到。


%===================================================================%



上述资料基本上完整地描述了原始问题,不过有几点内容附带说明一下:



1. 首先更正一个错误,我认为在“生成N ( 0.0128, 0.9596 ) 的高斯分布序列”的程序中,应该改为以下的代码:


%===================eight=====================================% 
 
 y=randn(1,2500); 
 
 y=y-mean(y); 
 
 y=y/std(y); 
 
 a=0.0128; 
 
 b=sqrt(0.9596); 
 
 y=a+b*y; 
 
 %==========================================================%


2. 上面资料最后部分隐含了一个出自zhyuer 版友的结论:


%==========================zhyuer===================================%


1)        rand产生的是[0,1]上的均匀分布的随机序列


2)        randn产生均值为0,方差为1的高斯随机序列,也就是白噪声序列;


%===================================================================%


也就是说,可以直接使用上面两个函数对原始信号添加噪声(例如y=x+rand(length(x),1)或者y=x+randn(length(x),1))



3. 事实上,无论是wgn还是awgn函数,实质都是由randn函数产生的噪声。即,wgn函数中调用了randn函数,而awgn函数中调用了wgn函数。下面就我熟悉的“向已知信号添加某个信噪比(SNR)的高斯白噪声”来说明一下,不过如果大家阅读过awgn的实现代码就不用看下去了,呵呵。从上述可知,这个任务可以使用awgn函数实现,具体命令是:awgn(x,snr,’measured’,'linear’),命令的作用是对原信号f(x)添加信噪比(比值)为SNR的噪声,在添加之前先估计信号f的强度。这里涉及三个问题:在awgn这个函数中,SNR是如何计算的?什么是信号的强度?awgn函数具体是如何添加噪声的?事实上,前两个问题是相关的,因为根据定义,SNR就是信号的强度除以噪声的强度,所以,首先来讲讲信号的强度。其实信号的强度指的就是信号的能量,在连续的情形就是对f(x)平方后求积分,而在离散的情形自然是求和代替积分了。在matlab中也是这样实现的,只不过多了一个规范化步骤罢了:


sigPower = sum(abs(sig( 
 
).^2)/length(sig( 
 
)


这就是信号的强度。至此,SNR的具体实现也不用多说了(注:

由于采用的是比值而非db,所以与下面“计算信噪比”所使用的方式不同,即没有求对数步骤)。

最后说说awgn函数具体是如何添加噪声的。事实上也很简单,在求出f的强度后,结合指定的信噪比,就可以求出需要添加的噪声的强度noisePower=sigPower/SNR。由于使用的是高斯白噪声即randn函数,而randn的结果是一个强度为1的随机序列(自己试试sum(randn(1000,1).^2)/1000就知道了,注意信号的长度不能太小)。于是,所要添加的噪声信号显然就是:sqrt(noisePower)*randn(n,1),其中n为信号长度。



4. 上面所说的都是具有分布特性(相关的)随机序列,如果需要添加不相关的随机序列,则可以使用jimin版友的方法:


%========================jimin=======================================% 
 
 for i=1:100 
 
 x(i)=randn(1); 
 
 end 
 
 %===================================================================%


即先产生噪声信号,后再与原信号叠加。




最后是另外的一些常见问题,整理如下:



1. Matlab中如何产生值为0,1的随机序列?【转bainhome版友】:round(rand(5))



2. Matlab中如何计算信噪比?下面的代码转自Happy教授:



%===========================Happy===================================% 
 
 function snr=SNR(I,In) 
 
 % 计算信号噪声比函数 
 
 % by Qulei 
 
 % I 
 
riginal 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. 随机产生1-n的索引排列:randperm函数