2013-10-03 134 views
-1

嘿傢伙我試圖找到我記錄的.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條記錄

任何想法???由於

+0

L = x(i *(N-1):i * N); – durasm

+0

它只是意味着你的輸入'x'沒有那麼多元素。檢查函數第一行中的'size(x)'並將其發佈到此處。 –

回答

1

這顯然是家庭作業,所以我不會那樣做你的工作你。但是你的代碼有很多錯誤。通過固定所有這些首先開始:

  • 使用更合適的函數名,homework123不是一個好名字來描述函數的功能。

  • 使用更合適的變量名。在這方面更多的標準是nfft而不是Nn_average而不是NumRec。我不關心你使用的確切的東西,但它應該準確地描述變量的作用。

  • 你的錯誤信息明確暗示你正試圖指數x在一些非法的方式。首先創建一個只打印正確索引(1..1024,1025..2048,...)的循環,並確保它遵循指令E.只有當它按預期工作時,纔會添加其餘的代碼。

  • 您使用三重嵌套的循環。您只需要一個for-loop或while循環來解決這個問題。

+0

謝謝你的建設性批評:)我自己想出來了,不過謝謝你的提示 – user2643355