matlab中噪声有关问题

原文链接：http://blog.sina.com.cn/s/blog_9315d81e0102vin3.html

以matlab中awgn函数为例说明：

    在matlab中无论是wgn还是awgn函数，实质都是由randn函数产生的噪声。即：wgn函数中调用了randn函数，而awgn函数中调用了wgn函数。

    根据awgn的实现代码可以知道“向已知信号添加某个信噪比（SNR）的高斯白噪声”，即：awgn(x,snr,’measured’,'linear’)，命令的作用是对原信号x添加信噪比（比值）为SNR的噪声，在添加之前先估计信号x的强度。

    这里涉及三个问题：在awgn这个函数中，SNR是如何计算的？什么是信号的强度？awgn函数具体是如何添加噪声的？

    事实上，前两个问题是相关的，因为根据定义，SNR就是信号的强度除以噪声的强度（或者信号功率与噪声功率之比），所以，首先来讲讲信号的强度。其实信号的强度指的就是信号的能量，在连续的情形就是对x平方后求积分，而在离散的情形自然是求和代替积分了。在matlab中也是这样实现的，只不过多了一个规范化步骤罢了：

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

这就是信号的强度，这里sig(: )为信号。

    至此，SNR的具体实现也不用多说了（注：由于采用的是比值而非db，所以与下面“计算信噪比”所使用的方式不同，即没有求对数步骤）。

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

自然要求的白噪声的方差也可以求出来了!

更新程序如下：

function snr=SNR(I,In)

% 计算信号噪声比函数

% 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

*******************************************************************************************************************

一个例子：

X = sqrt(2)*sin(0:pi/1000000:6*pi);         %产生正弦信号

Y = awgn(X,10,'measured');                  %加入信噪比为10db的噪声，加入前预估信号的功率（强度）

sigPower = sum(abs(X).^2)/length(X)         %求出信号功率

noisePower=sum(abs(Y-X).^2)/length(Y-X)     %求出噪声功率

SNR=10*log10(sigPower/noisePower)           %由信噪比定义求出信噪比，单位为db

*******************************************************************************************************************

关于wgn

WGN（m，n，p）产生功率为p dBW的m*n的高斯白噪声矩阵，其中p是以dbW为单位的输出强度。 

若要产生一个均值0，方差为0.0965 的高斯白噪声，不可直接用WGN（N，1，0.0965）产生，而应该如下：

1. N=1000;
    x=sqrt(0.0965)*randn(N,1);
    Px=(x.'*x)/N   % 验证，这里Px的求法与上面noisePower=sum(abs(Y-X).^2)/length(Y-X)的求法是一致的

2. N=1000;
    y=wgn(N,1,10*log10(0.0965));
    Py=(y.'*y)/N   % 验证

一点说明，对高斯白噪声，其方差和功率（单位为W）是一样的。因此，对方差，要做的只是将w变换成dbw，即dbw=10log（w）。

*******************************************************************************************************************

信噪比，英文名称叫做SNR或S/N（Signal Noise Ratio)，是指系统中信号与噪声的比例。信号指的是来自设备外部需要通过这台设备进行处理的电子信号，噪声是指经过该设备后产生的原信号中并不存在的无规则的额外信号（或信息），并且该种信号并不随原信号的变化而变化。

信噪比的计量单位是dB，其计算方法是10LOG(Ps/Pn)，其中Ps和Pn分别代表信号和噪声的有效功率，也可以换算成电压幅值的比率关系：20LOG(Vs/Vn)，Vs和Vn分别代表信号和噪声电压的“有效值”。信噪比应该越高越好。

原文网址：http://www.ilovematlab.cn/forum.php?mod=viewthread&tid=54155

在信号处理中经常需要把噪声叠加到信号上去，在叠加噪声时往往需要满足一定的信噪比，这样产生二个问题，其一噪声是否按指定的信噪比叠加，其二怎么样检验带噪信号中信噪比满足指定的信噪比。
在MATLAB中可以用randn产生均值为0方差为1的正态分布白噪声，但在任意长度下x=randn(1,N)，x不一定是均值为0方差为1（有些小小的偏差），这样对后续的计算会产生影响。在这里提供3个函数用于按一定的信噪比把噪声叠加到信号上去，同时可检验带噪信号中信噪比。
1，把白噪声叠加到信号上去：
function [Y,NOISE] = noisegen(X,SNR)
% noisegen add white Gaussian noise to a signal.
% [Y, NOISE] = NOISEGEN(X,SNR) adds white Gaussian NOISE to X.  The SNR is in dB.
NOISE=randn(size(X));
NOISE=NOISE-mean(NOISE);
signal_power = 1/length(X)*sum(X.*X);
noise_variance = signal_power / ( 10^(SNR/10) );
NOISE=sqrt(noise_variance)/std(NOISE)*NOISE;
Y=X+NOISE;
其中X是纯信号，SNR是要求的信噪比，Y是带噪信号，NOISE是叠加在信号上的噪声。

2，把指定的噪声叠加到信号上去
有标准噪声库NOISEX-92，其中带有白噪声、办公室噪声、工厂噪声、汽车噪声、坦克噪声等等，在信号处理中往往需要把库中的噪声叠加到信号中去，而噪声的采样频率与纯信号的采样频率往往不一致，需要采样频率的校准。
function [Y,NOISE] = add_noisem(X,filepath_name,SNR,fs)
% add_noisem add determinated noise to a signal.
% X is signal, and its sample frequency is fs;
% filepath_name is NOISE's path and name, and the SNR is signal to noise ratio in dB.
[wavin,fs1,nbits]=wavread(filepath_name);
if fs1~=fs
    wavin1=resample(wavin,fs,fs1);
end
nx=size(X,1);
NOISE=wavin1(1:nx);
NOISE=NOISE-mean(NOISE);
signal_power = 1/nx*sum(X.*X);
noise_variance = signal_power / ( 10^(SNR/10) );
NOISE=sqrt(noise_variance)/std(NOISE)*NOISE;
Y=X+NOISE;
其中X是纯信号，filepath_name是指定噪声文件（.wav）的路径和文件名，SNR是要求的信噪比，fs是信号X的采样频率，Y是带噪信号，NOISE是叠加在信号上的噪声。

3，检验带噪信号的信噪比
信噪比的定义为
           信号能量              (纯信号)^2
SNR=-----------------=--------------------------
           噪声能量        (带噪信号-纯信号)^2

function snr=SNR_singlech(I,In)
% 计算信噪比函数
% I :original signal
% In:noisy signal(ie. original signal + noise signal)
snr=0;
Ps=sum(sum((I-mean(mean(I))).^2));%signal power
Pn=sum(sum((I-In).^2));           %noise power
snr=10*log10(Ps/Pn);
其中I是纯信号，In是带噪信号，snr是信噪比

以下给出调用上函数的例子可作参考：
例一
clear all; clc; close all;
[filename,pathname]=uigetfile('*.wav','请选择语音文件：');
[X,fs]=wavread([pathname filename]);
[Y,NOISE] = noisegen(X,10);
subplot 311; plot(X);
subplot 312; plot(NOISE);
subplot 313; plot(Y);
mn=mean(NOISE)
snr=SNR_singlech(X,Y)

例二
clear all; clc; close all;
[filename,pathname]=uigetfile('*.wav','请选择语音文件：');
[filename1,pathname1]=uigetfile('*.wav','请选择噪声文件：');
filepath_name=[pathname1 filename1];
[X,fs]=wavread([pathname filename]);
[Y,NOISE] = add_noisem(X,filepath_name,10,fs);
subplot 311; plot(X);
subplot 312; plot(NOISE);
subplot 313; plot(Y);
mn=mean(NOISE)
snr=SNR_singlech(X,Y)