MATLAB的答案

关于matlab pWelch实现的问题

135次浏览(过去30天)
显示旧的评论
Neraj
Neraj 2012年2月20日
  • 0
编辑: 摩西2019年5月8日
你好,
我有两个关于matlab的pWelch实现的问题。如果我的问题被别人问了，我很抱歉，但我已经花了很长时间寻找我的问题的答案，我还没有偶然发现他们。所以我决定开始一个新的问题。
一些背景知识:我是用matlab来做原型的。现在我想要在c#中部署我的代码，所以我必须从matlab中翻译一切。在我采取这一步之前，我想我应该首先将matlab的pWelch分解成它的“子步骤”，并确保我得到与原始函数相同的结果。我很接近，但没有得到相同的结果。根据pWelch的文档，mathworks的工作人员遵循以下步骤:
1.将输入信号向量x按窗口和重叠(或其默认值)划分为k个重叠段。如果窗口大小大于FFT点(NFFT)的数量，则将信号分成NFFT长度的段，然后将最后一个段填充为0。
2.指定的(或默认的)窗口应用于x的每个段。(在将窗口应用于每个段之前不做预处理。)
3.一个非FFT点FFT应用于有窗数据。
4.计算每个加窗段的(修改)周期图。
5.对修正的周期图集进行平均，形成谱估计S(jω)。
6.得到的频谱估计被缩放以计算功率谱密度为S(jω)/F，其中F是-2π当你不提供采样频率-fs当你提供采样频率
所以我只是遵循他们的每一步，一个接一个，来重现结果。然而，我偏离了一些“标准化因素”。我的结果都乘以这个因数。另一个区别是:0Hz, 1Hz和256Hz组件是错误的。其他都是正确的。这真的很令人困惑。
我深入挖掘了一下，发现了这篇文章: < //www.tatmou.com金宝app/support金宝搏官方网站/solutions/en/data/1-17M0Q/index.html?product=SG&solution=1-17M0Q >
我并没有真正掌握整个“有偏vs无偏”的概念，但我试着乘以范数(w)^2/和(w)^2，它仍然没有解决我的问题。我注意到他们引用了一本书，第419页 http://i.imgur.com/7mvak.png
所以我看了一下那本书，我发现他们确实有一个不同的标准化因子。我也试过了，还是没找到。
下面是我编写的比较代码。我的数据采样在512Hz。我正在做一个窗口长度512点的pWelch, FFT长度512点，有50%的重叠。使用标准汉明窗。我用了汉德尔作为假数据，这样你们就能更容易地运行代码。
负载汉德尔；
数据= y (1:2048) ';
做一个汉明窗
w =汉明(512);
%为PSD结果预分配空间
mypsd = 0 (512 7);
%步骤1:遍历数据，每次512个点，256个点重叠
为i = 0:6
步骤2:应用一个汉明窗口
临时数据=(1 + 256 *我:我* 256 512 +)”。* w;
%步骤3:计算FFT，取前半部分
temp = fft(临时);
第四步:计算周期图的绝对值的平方
temp = abs(临时)^ 2;
将结果保存在存储变量中
mypsd(:,我+ 1)= temp;
结束
第五步:对所有周期图取平均值
mypsd_v1 =意味着(mypsd, 2);
%步骤6:除以采样率(512Hz为我的数据)
mypsd_v1 = mypsd_v1/512;
扔掉mypsd的后半部分
mypsd_v1 = mypsd_v1 (1:257);
%的matlab等效，使用pWelch函数
[p f] = pwelch(数据,512年,.5,512,512);
%现在计算乘数据的因子，所以我的方法和
用matlab的方法进行匹配
norm_factor = mypsd_v1 (115) / p (115);
mypsd_v1 = mypsd_v1 / norm_factor;
%绘制结果
情节(0:长度(p) 1、日志(p),“b”),持有在、情节(0:长度(p) 1、日志(mypsd_v1),“r”）;
标题(sprintf (“归一化因素:% f”norm_factor));
包含(的频率(赫兹)）;
ylabel (“日志实力”）;
%版本2:应用U因子。L =窗口长度，K =截面的#
U =(1/512) *总和(w。^ 2);
K = 7;
L = 512;
mypsd_v2 = (1 / (K * L * U)) * (mypsd, 2)总和;
mypsd_v2 = mypsd_v2 (1:257);
norm_factor = mypsd_v2 (115) / p (115);
mypsd_v2 = mypsd_v2 / norm_factor;
%绘制结果
人物,情节(0:长度(p) 1、日志(p),“b”),持有在、情节(0:长度(p) 1、日志(mypsd_v2),“r”）;
标题(sprintf (“归一化因素:% f”norm_factor));
包含(的频率(赫兹)）;
ylabel (“日志实力”）;
任何帮助/指示将非常感谢!
谢谢,Neraj

登录回答这个问题。

答案(4)

摩西
摩西 2016年9月20日
编辑:摩西 2019年5月8日
嗨Neraj,
你的代码能正常工作吗?子步骤是正确的，但代码的某些部分似乎出错了，特别是规范化因子。我对您的代码做了一些更改，所以现在显示了相同的结果。
关于为什么使用这个归一化因子，请参考本文档中的公式#24。这是一篇非常有用和详细的关于谱密度估计和窗口的论文。 http://holometer.fnal.gov/GH_FFT.pdf
2来自于忽略多余的负频率。除以fs是因为想要得到功率谱 密度 表达了在V ^ 2 /赫兹。窗口的平方和来自于两个光谱相互相乘的结果(Pxx = Px⋅Px*，其中⋅代表点积，∗代表复共轭)。希望这也能帮助未来的观众了解这个答案。
负载汉德尔；
数据= y (1:2048) ';
做一个汉明窗
fs = 512;
w =汉明(512);
%为PSD结果预分配空间
Mypsd = 0 (512, 7);
%步骤1:遍历数据，每次512个点，256个点重叠
为i = 0:6
步骤2:应用一个汉明窗口
临时数据=(1 + 256 *我:我* 256 512 +)”。* w;
%步骤3:计算FFT，取前半部分
temp = fft(临时);
第四步:计算周期图的绝对值的平方
temp = abs(临时)^ 2;
将结果保存在存储变量中
Mypsd (:， i+1) = temp;
结束
第五步:对所有周期图取平均值
mypsd_v1 =意味着(mypsd, 2);
扔掉mypsd的后半部分
mypsd_v1 = mypsd_v1 (1:257);
%规格化因素
mypsd_v1 = mypsd_v1 / (fs *总和(w ^ 2));
%忽略DC和Nyquist值
Mypsd_v1 (2:end-1) = Mypsd_v1 (2:end-1) * 2;
%的matlab等效，使用pWelch函数
[pxx, f] = pwelch(data, w， []， 512, fs);
图;
次要情节(2,1,1);
情节(f, mypsd_v1);
次要情节(2,1,2);
情节(f, pxx);
2的评论
显示 1年长的评论
沃尔特·罗伯森
沃尔特·罗伯森 2018年4月30日
四个512的完整窗口使得长度为2048。在半窗口偏移量处加上三个窗口，总共7个。最后一个以0为基础的索引是6。

登录并发表评论。

韦恩王
韦恩王 2012年2月20日
你好，Neraj，你的问题和代码很长，所以我只是给你一个简单的例子，使用两个窗口，没有重叠的随机向量。我先把随机数生成器设为它的默认值，这样我就可以给出最终结果，但这并不重要。您将看到以下内容与MATLAB pwelch()一致。
rng默认的；
x = randn (1024 1);
Pxx = pwelch (x, 512, 0512);
水稻播种期及秧龄= 0 (512 2);
为k = 0:1
x1 = x (1 + k * 512: (k + 1) * 512);
xw = x1。*汉明(512);
水稻播种期及秧龄(:,k + 1) = abs (fft (xw)) ^ 2. /规范(汉明(512),2)^ 2;
结束
水稻播种期及秧龄=意味着(水稻播种期及秧龄、2);
水稻播种期及秧龄= xper. /(2 *π);
水稻播种期及秧龄=水稻播种期及秧龄(1:长度(水稻播种期及秧龄)/ 2 + 1);
水稻播种期及秧龄(2:end-1) = 2 *水稻播种期及秧龄(2:end-1);
马克斯(abs (Pxx-xper))
可以看到绝对值的最大差值只有10^(-16)的数量级。
希望这有助于
2的评论
显示 1年长的评论
彼得
彼得 2018年11月7日
嗨，韦恩，我试图推广以上代码的情况下，段长度和nfft是不同的成功的情况下，nfft大于段长度。参见下面的情形2。但是，nfft小于段长度的情况(情况1)让我很困惑。Matlab 2017b帮助说:“如果nfft小于段长度，段被包装使用datawrap使长度等于nfft。”但我不能理解这是什么意思。你能照点光吗?(可能通过一段代码的方式，在问号所在的地方)
清晰的所有
casenr = 2
元= 1024;
开关casenr
情况下1比nfft长%的片段，2个片段
segment_length = nt / 2;
nfft = 400;
情况下2%段小于nfft, 2段
segment_length = 400;
nfft = 512;
结束
nzeropad = nfft-segment_length;
x = randn (nt, 1);
[Pxx w] = pwelch (x, segment_length 0 nfft);
水稻播种期及秧龄= 0 (nfft, 2);
为k = 0:1
%将段
i1 = 1 + k * segment_length;
i2 = (k + 1) * segment_length;
x1 = x (i1: i2);
%窗口+零填充:
如果nfft > = segment_length
windowlength = segment_length;
xw = [x1。*汉明(windowlength); 0 (nzeropad 1)];
elseifnfft < segment_length
%?????????
结束
水稻播种期及秧龄(:,k + 1) = abs (fft (xw)) ^ 2. /规范(汉明(windowlength), 2) ^ 2;
结束
水稻播种期及秧龄=意味着(水稻播种期及秧龄、2);
水稻播种期及秧龄= xper. /(2 *π);
水稻播种期及秧龄=水稻播种期及秧龄(1:长度(水稻播种期及秧龄)/ 2 + 1);
水稻播种期及秧龄(2:end-1) = 2 *水稻播种期及秧龄(2:end-1);
马克斯(abs (Pxx-xper))

登录并发表评论。

韦恩王
韦恩王 2012年2月22日
你复制粘贴了这段代码?
x = randn (1024 1);
Pxx = pwelch (x, 512, 0512);
水稻播种期及秧龄= 0 (512 2);
为k = 0:1
x1 = x (1 + k * 512: (k + 1) * 512);
xw = x1。*汉明(512);
水稻播种期及秧龄(:,k + 1) = abs (fft (xw)) ^ 2. /规范(汉明(512),2)^ 2;
结束
水稻播种期及秧龄=意味着(水稻播种期及秧龄、2);
水稻播种期及秧龄= xper. /(2 *π);
水稻播种期及秧龄=水稻播种期及秧龄(1:长度(水稻播种期及秧龄)/ 2 + 1);
水稻播种期及秧龄(2:end-1) = 2 *水稻播种期及秧龄(2:end-1);
马克斯(abs (Pxx-xper))
首先清理你工作空间中的所有东西，然后尝试一下。
我完全同意。你用的是什么版本的MATLAB ?
4评论
显示 3年长的评论
Celso Coslop Barbante
Celso Coslop Barbante 2012年7月26日
还有协议。这段代码很好地理解了如何用FFT模拟pwelch PSD估计器。我试图实现一个类似pwelch的功能在硬件使用FPGA和FFT核心。谢谢。

登录并发表评论。

Invizible灵魂
Invizible灵魂 2012年11月12日
你好,
为什么在下面这行代码中
水稻播种期及秧龄= xper. /(2 *π);
这里有一个除法 2 *π ？
诚恳地。

类别

标签

社区寻宝

在MATLAB中央找到宝藏，发现社区如何可以帮助你!

开始狩猎!