加速度信号中加10%的噪声

secondye

上一个帖子问，怎么在Matlab中实现对桁架模型施加白噪声激励，得出加速度响应，然后对响应信号处理，得到结构模态参数。
感谢VibrationMaster提供的帮助。现在问题解决了。

现在又遇到个小问题，就是为了验证模态参数识别程序的有效性，打算在加速度响应信号里添加5%，10%，20%的噪声，看到文章上有人这样做的，我也想验证一下我的程序在噪声环境下是否能准确识别出参数。
现在就是不知道怎么样在已有响应信号中添加5%，10%，20%的噪声，请求一下各位。

感谢

zhouyang664

添加噪声好像有专门的函数吧！
网上搜索一下，或等待高人！

Happy99

help awgn?
个人以为那只是说明一种大小比较而已
ref: 关于信噪比与向已知信号添加噪声的整理及个人心得http://forum.vibunion.com/thread-31086-1-1.html

VibrationMaster

1. 假定x(1),x(2),....,x(n)是计算得到的加速度.
2. 计算出它的方差 sigma
3. 用 sigma*sqrt(10%)*randn(1,n) 或 sigma*sqrt(10%)*randn(n,1)产生n个随机数
4.将n个随机数分别加到x(1),x(2),....,x(n)-->完成任务

secondye

回复 4 # VibrationMaster 的帖子

再次感谢 VibrationMaster的解答。
刚刚试了一下，加噪声前后做FFT分析的结果都很好，峰值都很明显，按道理应该加噪声后的FFT结果不是那么好吧？
不知道代码跟方法有没有写错，贴上了再请教：
% gAcce是原始响应信号，gAcce_noise是添加噪声的响应信号
% Nodf 是响应通道数目,npts是数据点数目
for i=1:Ndof
    sigma=std(gAcce,0,1);                         %% 求每一列信号的均方差
    wn(:,i)=randn(npts,1);                         %% 产生均值为0，方差为1的高斯随机序列，也就是白噪声序列
    noise(:,i)=sigma(i)*sqrt(0.1)*wn(:,i);    %% 按照VibrationMaster的方法
    gAcce_noise(:,i)=gAcce(:,i)+noise(:,i);  %% 信号叠加
end

              原始信号做FFT            
  
       加10%噪声的信号做FFT

VibrationMaster

FFT本质相当于一种平均，它相当于达到强化周期信号的谱峰效果，而白噪声能量则是分布在所有的频率带上。--这也就是为什么谱看起来比时程曲线看起来更漂亮的原因，也是用的比较广的原因--信号/噪声分离的比较直观

ChaChing

VibrationMaster 发表于 2010-11-6 20:06

                               
登录/注册后可看大图

...用 sigma*sqrt(10%)*randn ...

好奇为何需取sqrt ?

secondye

回复 6 # VibrationMaster 的帖子

哦，正是，添加白噪声后的频谱图，明显看出底下的应该就是白噪声的能量分布吧？这在添加白噪声前的FFT图上是没有的。

我想问的是，一般的实际信号采集过程中，有环境噪声的影响，对信号进行分析时，一般会出现一些虚假的模态，对识别真实模态的精度造成影响。
我想模拟这个信号采集过程。比较无噪声和有噪声情况下，我的程序识别模态参数的精度。前面我做的添加白噪声的过程能不能当作是添加环境噪声？
如果是这样的一个过程，如您所说的，白噪声的能量分布在所有的频率带上，对FFT分析精度没有影响？那对时域法如ERA，SSI的识别精度呢？

VibrationMaster

回复 7 # ChaChing 的帖子

噪声一般按能量计算，换算成幅度应开个平方

VibrationMaster

回复 8 # secondye 的帖子

1。所有的估计算法都会受到噪声的影响，但是国内没有人做这个系统的比较。你感兴趣可以做个硕士毕业论文。
2。按照统计的观点，应该生成N组随机噪声，然后统计这N组的模态参数的均值和方差，用这组均值和方差来评估算法
3。上述的N应该大于1000

ChaChing

回复 9 # VibrationMaster 的帖子

可能我误解了:@L
所谓的方差指的是variance(var)而不是standard deviation(std)!?

ChaChing

VibrationMaster 发表于 2010-11-7 14:37

                               
登录/注册后可看大图

回复 8 # secondye 的帖子

1。所有的估计算法都会受到噪声的影响，但是国内没有人做这个系统的比较。你感兴 ...

所有的估计算法都会受到噪声的影响应该容易理解
但比较好奇的是, 国内没有人做过系统的比较!? 没教授有兴趣?

VibrationMaster

确实是这样，国际上各人都说自己的方法好，但是系统地用噪声干扰来比较的文献我还没有读到，包括英文。
在发表论文界，这是吃力不讨好的活，尤其在振动界。我现在的主要方向不做这个了，想找个学生把所有的算法都读懂并实现，但是现在学生太懒，太ben，而且也没有这方面的基金支撑

secondye

回复 13 # VibrationMaster 的帖子

哦，刚刚在Matlab中文论坛看到一个好帖，关于叠加噪声的，转过来了‘
原帖地址：http://www.ilovematlab.cn/thread-54155-1-1.html

Matlab信号上叠加噪声和信噪比的计算
在信号处理中经常需要把噪声叠加到信号上去，在叠加噪声时往往需要满足一定的信噪比，这样产生二个问题，其一噪声是否按指定的信噪比叠加，其二怎么样检验带噪信号中信噪比满足指定的信噪比。
在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)

Happy99

回复 14 # secondye 的帖子

3F连接老八的帖子不是有吗!?:@)