嘿傢伙我試圖找到我記錄的.wav信號的功率譜密度,這實際上是一個沉浸在噪聲中的正弦波。我寫的函數應該記錄全部1024個點的記錄,並通過查找每個記錄的Gxx,然後添加它們並將它們除以在下面的算法中更好地解釋的記錄數來找到信號的Gxx :???索引超過矩陣尺寸PSD Proplem
a。瀏覽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
我不確定h解決它,我已經嘗試了許多不同的方法,但仍然得到這個錯誤。我對Matlab不太熟悉,但是我真正想爲第19行做的所有事情就像:
x(1:1024),x(1025:2048),x(2049:3072),x(3072 :4096)...等等到100條記錄
任何想法???由於
L = x(i *(N-1):i * N); – durasm
它只是意味着你的輸入'x'沒有那麼多元素。檢查函數第一行中的'size(x)'並將其發佈到此處。 –