？？？索引超出矩阵维度 PSD 问题

Question

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

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

b.对此记录执行正常的 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 条记录

有什么想法吗？？？谢谢

Answer 1

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

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

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

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

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