Matlab高斯白噪声通过互相关函数求功率谱密度 |
您所在的位置:网站首页 › 求平均功率的两种方法 › Matlab高斯白噪声通过互相关函数求功率谱密度 |
参考文章:百度文库文章 但是该文章有很多部分,不知是matlab更新了还是咋地,有错误,踩了很多坑,仅做参考。 产生高斯白噪声的方法在“信号处理”专栏中有写,其时域和频域图如下: 用matlab中求自(互)相关的xcorr函数,参考帮助文档 [r,lags]=xcorr(z,'biased');其中z是上述高斯白噪声,r是自相关函数,lags是时间偏移量(索引),尤其注意’biased’参数,这是调试了半天才发现的问题。 帮助文档里描述如下: 相当于算出来的r除以了样本数,这是离散自相关函数计算步骤。 如果不加这个参数,自相关函数值会变得很大。 算出自相关函数之后,根据维纳辛钦定理,信号的功率谱密度和信号的自相关函数是一个傅里叶变换对。使用fft函数求自相关函数的傅里叶变换: pdv=fft(r); pdv=abs(fftshift(pdv))./length(pdv);fftshift是频谱矫正函数,将fft后的结果以中心为分界,两边分别做镜像,得出来的结果再除以样本数,即双边频谱图。 而这个频谱图就是原信号的功率谱密度。 plot((0:length(pdv)-1)*L/length(pdv)-L/2,pdv)注意这里横坐标是频率,需要做映射处理,L为样本数量 高斯白噪声其实看不出来啥,下面用一个没加噪声的纯信号来展示: y=12*cos((2*pi)*100.*t)+15*cos((2*pi)*150.*t)+18*cos((2*pi)*210.*t)+10;上面的功率谱是由双边频谱的平方(除直流)计算出来的,下面的是用维纳辛钦定理计算出来的,可以看到功率谱几乎相同。 信号的自相关函数如下: 编者水平有限,很多知识点都是似懂非懂,如有错误欢迎指出! 附代码: fs=1000;%采样频率hz T_N=1.5;%总时间s t=1/fs:1/fs:T_N;%时间向量 L=T_N*fs;%样本长度 y=12*cos((2*pi)*100.*t)+15*cos((2*pi)*150.*t)+18*cos((2*pi)*210.*t)+10;%信号 subplot(2,1,1); plot(t,y); xlabel("时间/s") ylabel("幅度/v") title("时域") fft_y=fft(y);%快速傅里叶变换 P = abs(fft_y/L);%取幅频特性,除以L P = P(1:L/2+1);%截取前半段 P(2:end-1)=2*P(2:end-1);%单侧频谱非直流分量记得乘以2 f = fs*(0:(L/2))/L;%频率,最多到一半(奈奎斯特采样定理) subplot(2,1,2); plot(f,P); xlabel("频率/Hz") ylabel("幅度/v") title("单边频谱") figure(4) [r,lags]=xcorr(y,'biased');%得到自相关函数的幅度和偏移量 subplot(2,1,1) plot(lags,r) xlabel("时间偏移/s") ylabel("相关程度") title("(自相关函数)") fft_y=fftshift(fft_y);%频谱矫正 powerpu=abs(fft_y./L).^2;%双边功率谱 % subplot(2,1,1) % plot((1:size(powerpu,2))-751,powerpu) % xlabel("频率/Hz") % ylabel("功率/W") % title("功率谱") pdv=fft(r,size(r,2));%对自相关函数快速傅里叶变换 pdv=abs(fftshift(pdv));%频谱矫正,让正半轴部分和负半轴部分的图像分别关于各自的中心对称,得到双边谱 subplot(2,1,2) plot((1:length(pdv))*L/length(pdv)-L/2,pdv./length(powerpu)) xlabel("频率/Hz") ylabel("功率谱W/Hz") title("(功率谱)") % z1=0.1*randn(1,201);%产生方差N(0,0.12)高斯白噪声 % [r1,lags]=xcorr(z1,'unbiased');%自相关函数的估计 % plot(lags,r1); % f1=fft(r1); % f2=fftshift(f1);%频谱校正 % l1=((0:length(f2)-1)*200/length(f2)-100);%功率谱密度x轴 % y4=abs(f2); % figure(2) % plot(l1,y4); |
今日新闻 |
点击排行 |
|
推荐新闻 |
图片新闻 |
|
专题文章 |
CopyRight 2018-2019 实验室设备网 版权所有 win10的实时保护怎么永久关闭 |