-1

嘿伙计们，我正在尝试找到我记录的 .wav 信号的功率谱密度，该信号本质上是沉浸在噪声中的正弦波。我编写的函数应该记录所有 1024 个点的长度，并使用它通过查找每条记录的 Gxx 来找到信号的 Gxx，然后将它们相加并将它们除以在下面的算法中更好地解释的记录数：

一个。遍历 wav 文件并提取第一条记录长度（例如 1 到 1024 点）。（请注意，记录长度是您的新“N”，因此频率间隔会根据此变化，而不是 wav 文件的总长度）。

湾。对此记录执行正常的 PSD 功能。

C。存储此向量。

d。提取 wav 文件中接下来的 1024 个点（例如 1025:2048）并对此记录执行 PSD。

e. 将此添加到先前存储的记录中，并继续执行步骤 c 到 e，直到到达 wav 文件的末尾或所需的记录总数。（记住总记录*记录长度必须小于 wavfile 的总长度！）

F。将 PSD 除以平均数（或记录数）。这是你的平均 PSD

我创建的函数如下：

%Function to plot PSD
function[f1, GxxAv] = HW3_A_Fn_811003472_RCT(x,fs,NumRec)
Gxx = 0;
GxxAv = 0;

N = 1024;
df = fs/N;
f1 = 0:df:fs/2;
dt = 1/fs;
T = N*dt;

q = 0;
e = 1;

for i = 1:NumRec;
    for r = (1+q):(N*e);

        L = x(1+q:N*e);
        M = length(L);
        Xm = fft(L).*dt;
        aXm = abs(Xm);

        Gxx(1)=(1/T).*(aXm(1).*aXm(1));
        for k = 2:(M/2); 
            Gxx(k) = (2/T) *(aXm(k).*(aXm(k)));
            %Gxx = Gxx + Gxx1(k);
        end
        Gxx((M/2)+1)= (1/T)*(aXm((M/2)+1)).*(aXm((M/2)+1));
        q = q+1024;
        e = e+1;
        %Gxx = Gxx + Gxx1((M/2)+1);
    end
    GxxAv = GxxAv + Gxx;
    %Gxx = Gxx + Gxx1;
end
GxxAv = GxxAv/NumRec;

而我用来调用这个函数的代码如下：

[x,fs] = wavread('F:\Final\sem1Y3\Acoustics\sinenoise5s.wav');

[f1,GxxAv] = HW3_A_Fn_811003472_RCT(x,fs,100); %where 100 is the number of records to generated

plot(f1,GxxAv)

xlabel ('Frequency / Hz', 'fontsize', 18)
ylabel ('Amplitude Squared per Frequency / WU^2/Hz', 'fontsize', 18)
title ('Plot of the single sided PSD, using Averaging', 'fontsize', 18)
grid on

尝试绘制此图时，观察到以下错误：

??? Index exceeds matrix dimensions.

Error in ==> HW3_A_Fn_811003472_RCT at 19
        L = x(1+q:N*e);

Error in ==> HW3_A_3_811003472_RCT at 3
[f1,GxxAv] = HW3_A_Fn_811003472_RCT(x,fs,100); %where 100 is the number of records to generated

我不知道如何修复它，我尝试了许多不同的方法，但我仍然收到此错误。我对 Matlab 不太熟悉，但我真正想做的第 19 行就是：

x(1:1024), x(1025:2048), x(2049:3072), x(3072:4096)...等到 100 条记录

有任何想法吗？？？谢谢

4

1 回答 1

1

这显然是家庭作业，所以我不会为你做你的工作。但是您的代码有很多问题。首先修复所有这些：

  • 使用更合适的函数名称，homework123并不是描述函数作用的好名称。

  • 使用更合适的变量名。在这种情况下，更标准的做法是nfft代替N和n_average代替NumRec. 我不关心您使用的确切内容，但它应该准确描述变量的作用。

  • 您的错误消息清楚地暗示您正试图以x某种非法方式建立索引。从创建一个仅打印正确索引（1..1024、1025..2048，...）的循环开始，并确保它遵循您的指令 E。只有当它按预期工作时才添加其余代码。

  • 您使用三重嵌套的 for 循环。您只需要一个 for 循环或 while 循环即可解决此问题。

于 2013-10-03T19:19:49.620 回答