高分!急!用matlab分析功率谱密度,采样频率的设定.高手进!

Posted 2023-05-04

部分程序如下.解释一下,x_sample是我要分析的信号序列,sample rate是10Hz.我要分析该序列在[0.04 0.15]Hz区间和[0.15 0.4]区间的psd,并计算平均功率,考察它们的比(即r)随时间的变化情况.

我比较糊涂的是,psd设定的参数:
1.fs是采样频率,这里采样频率为什么要设成输入信号的频率?按照奈奎斯特定律,不是应该是它的两倍以上才可以包含全部频谱吗?
2.nfft是计算出来的点数,它的大小关系到后面得到的频域信号的精度,对吗?
3.单边oneside和双边除了幅度的区别还有什么别的?
4.单边的情况下y是4096(2^12)个点,那边对应的是0,那边是pi?

我最关心第一个问题,怎样才能确保我得到的频谱是正确的,否则后面计算的平均功率都扯淡了.如能给予解答,无限感激!

w=hann(m);
k=1;
Fs=10;
h = spectrum.periodogram('Hann');
m=100;
j=8;
for i=1:j:Npp-m
mx=(x_sample(i:m+i-1)-mean(x_sample(i:m+i-1)));
p=mx.*w';
hopts = psdopts(h);
set(hopts,'NFFT',2^13,'Fs',Fs,'SpectrumType','onesided');
Hpsd = psd(h,p,hopts);
y(:,k)=Hpsd.Data;
g1(k) = avgpower(Hpsd,[.04 .15]);
g2(k) = avgpower(Hpsd,[.15 .4]);
r(k)=g1(k)/g2(k);
time_r(k)=0.1*j*(k-1)+m/2*0.1;
k=k+1;
end;

楼主你好！

不知我的答案你是否满意：

（1）做波形显示以及fft变换，程序如下:

[y,fs]=wavread('E:\MATLAB6p5\work\3.wav');%读出信号，采样率。
y=y(:,1);%取单声道。
sigLength=length(y);
Y = fft(y,sigLength);
Pyy = Y.* conj(Y) / sigLength;
halflength=floor(sigLength/2);
f=Fs*(0:halflength)/sigLength;
figure;plot(f,Pyy(1:halflength+1));xlabel('Frequency(Hz)'); %画频域波形

t=(0:sigLength-1)/Fs;
figure;plot(t,y);xlabel('Time(s)'); %画时域波形

(2)关于滤波
声音频率主要集中在0~1KHZ，我想虑掉500hz以下的频率，因此采用一个高通滤波器
这里我使用了一个10阶butterworth高通滤波器，边带是500hz，但是这不能直接用，因为声音文件的采样率是44k，500hz相对于44k来说太小了。所以我得先把我的声音欠采样，然后再滤波。程序如下:
[k,Fs,bits]=wavread('E:\MATLAB6p5\work\3.wav');
k=k(:,1);
y_temp=k(1:90000);

dfactor=3;
y=decimate(y_temp,dfactor);
[b,a] = butter(10,500/(Fs/(dfactor*2)),'high');
y=filter(b,a,y);
y=interp(y,dfactor);

sigLength=length(y);
Y = fft(y,sigLength);
Pyy = Y.* conj(Y) / sigLength;
halflength=floor(sigLength/2);
f=Fs*(0:halflength)/sigLength;
figure;plot(f,Pyy(1:halflength+1));xlabel('Frequency(Hz)'); %滤波后的频域波形

t=(0:sigLength-1)/Fs;
figure;plot(t,y,t,y_temp);xlabel('Time(s)'); %滤波后的时域波形

wavwrite(y,Fs,'3_'); %保存处理后的声音文件，文件名为”3_”.

（3）
Matlab代码示例：

clear;
Fs=1000; %采样频率
n=0:1/Fs:1;

%产生含有噪声的序列
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));

window=boxcar(length(xn)); %矩形窗
nfft=1024;
[Pxx,f]=periodogram(xn,window,nfft,Fs); %直接法
plot(f,10*log10(Pxx));

间接法：
间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。

Matlab代码示例：

clear;
Fs=1000; %采样频率
n=0:1/Fs:1;

%产生含有噪声的序列
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));

nfft=1024;
cxn=xcorr(xn,'unbiased'); %计算序列的自相关函数
CXk=fft(cxn,nfft);
Pxx=abs(CXk);
index=0:round(nfft/2-1);
k=index*Fs/nfft;
plot_Pxx=10*log10(Pxx(index+1));
plot(k,plot_Pxx);

改进的直接法：
对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。

1. Bartlett法

Bartlett平均周期图的方法是将N点的有限长序列x(n)分段求周期图再平均。

Matlab代码示例：

clear;
Fs=1000;
n=0:1/Fs:1;
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
nfft=1024;
window=boxcar(length(n)); %矩形窗
noverlap=0; %数据无重叠
p=0.9; %置信概率

[Pxx,Pxxc]=psd(xn,nfft,Fs,window,noverlap,p);

index=0:round(nfft/2-1);
k=index*Fs/nfft;
plot_Pxx=10*log10(Pxx(index+1));
plot_Pxxc=10*log10(Pxxc(index+1));
figure(1)
plot(k,plot_Pxx);

pause;

figure(2)
plot(k,[plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc]);

2. Welch法

Welch法对Bartlett法进行了两方面的修正，一是选择适当的窗函数w(n)，并再周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时，可使各段之间有重叠，这样会使方差减小。

Matlab代码示例：

clear;
Fs=1000;
n=0:1/Fs:1;
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
nfft=1024;
window=boxcar(100); %矩形窗
window1=hamming(100); %海明窗
window2=blackman(100); %blackman窗
noverlap=20; %数据无重叠
range='half'; %频率间隔为[0 Fs/2]，只计算一半的频率

[Pxx,f]=pwelch(xn,window,noverlap,nfft,Fs,range);
[Pxx1,f]=pwelch(xn,window1,noverlap,nfft,Fs,range);
[Pxx2,f]=pwelch(xn,window2,noverlap,nfft,Fs,range);

plot_Pxx=10*log10(Pxx);
plot_Pxx1=10*log10(Pxx1);
plot_Pxx2=10*log10(Pxx2);

figure(1)
plot(f,plot_Pxx);

pause;

figure(2)
plot(f,plot_Pxx1);

pause;

figure(3)
plot(f,plot_Pxx2);

希望我的答案能够帮得到你！

谢谢！ 参考技术A 1.fs是采样频率,这里采样频率为什么要设成输入信号的频率?按照奈奎斯特定律,不是应该是它的两倍以上才可以包含全部频谱吗?
fs就是差值之后的序列的sample rate. 只要fs大于 sample的2倍以上频率就可以。
2.nfft是计算出来的点数,它的大小关系到后面得到的频域信号的精度,对吗?
是的。
3.单边oneside和双边除了幅度的区别还有什么别的?
没有别的了。
4.单边的情况下y是4096(2^12)个点,那边对应的是0,那边是pi?
左边是0右边是pi.本回答被提问者采纳 参考技术B 我是个路过挣分的 参考技术C 这个...以前做过...但是程序\题目\书都没带回来....
我的建议是你去搜索下"功率谱密度 matlab"
如果我没记错的话有不少资料的,这题目不难的
至于你的几个问题,再好好看看书吧,都是概念的东西

Matlab功率谱密度（PSD）实现方程

【中文标题】Matlab功率谱密度（PSD）实现方程

【英文标题】：Matlab power spectrum density(PSD) implement equation

【发布时间】：2017-09-06 01:13:09

【问题描述】：

伙计们，我正在学习如何使用 matlab 计算信号的 PSD。我知道函数Periodogram() 效果很好，但我想直接使用fft 方法。当我想详细学习时，我发现这个网站很有帮助：

https://www.mathworks.com/help/signal/ug/power-spectral-density-estimates-using-fft.html

但是，我在第一个示例中遇到了问题。我对这个说法有点困惑

psdx = (1/(Fs*N)) * abs(xdft).^2;

不知道为什么要用Fs和N作为参数。

【参考方案1】：

(1/(Fs*N)) 是一个归一化因子，以确保数值算法符合Parseval's Theorem。

【讨论】：

现在我明白为什么我们需要1/N。但是1/Fs 来自哪里？

@SihaoSun 1/Fs 从“power in bin”转换为“power per bin”，因此得到的psdx 是每个频率间隔的功率，而不是该频率附近的功率。

我完全理解你的意思。还有一个关于bin的小问题，bin是X-Y坐标轴上的一条线还是一个区域？

@SihaoSun 在此上下文中的“bin”是频率区间。严格来说，它不对应 X-Y 坐标平面中的任何东西，但我想如果你愿意的话，我想你可以把它想象成 X 轴上的一条线。