推 cawaiimi:推:) 10/24 00:35
※ 編輯: knuckles 來自: 125.230.31.100 (10/25 20:22)
※ 編輯: knuckles 來自: 125.230.227.152 (10/29 13:30)
我是有用fft做的
可能有人想知道怎麼做，所以來貼一下
首先STFT的數學式：
n+Q
X(m*df,n*dt) = Σ x(p*dt)exp(-j2π(m*df)(n*dt))dt
p=n-Q
也就是把陣列x的某一段(n-Q~n+Q共2*Q+1點)的值做FT後
存進矩陣X的第n行
直接寫成迴圈的方法：
X=zeros(length(f),length(t));
for n=1:length(t)
for m=1:length(f)
for p=n-Q:n+Q
X(m,n)=X(m,n)+x(p)*exp(-j*2*pi*(m-1)*df*(p-1)*dt)*dt;
end
end
end
注意陣列x與矩陣X的index一定要是大於等於1的整數
所以放進m、n、p就可以，不用乘dt、df
而後面的exp(...)項是數學式，所以index要乘dt、df
又矩陣X(m=1,n=1)時是用來代表數學式X(0,0)時的值
所以後面的數學式exp()項中出現的index要減1
不過這會產生一個問題，在n在開頭或結尾時
p=n-Q:n+Q會小於等於0或超過陣列x的長度
所以要把這種情況排除掉，改成下面這樣：
X=zeros(length(f),length(t));
for n=Q+1:length(t)-Q
for m=1:length(f)
for p=n-Q:n+Q
X(m,n)=X(m,n) + x(p)*exp(-j*2*pi*(m-1)*df*(p-1)*dt)*dt;
end
end
end
p就會落在1~length(t)之間了
接著改成使用fft的方式，fft的數學式：
N-1
X[m] = Σ x[k]*exp(-j2π*m*k/N), 1 <= m <= N.
k=0
程式的寫法為X=fft(x)，寫成迴圈就是：
X=zeros(1,N);
for m=1:N
for k=1:N
X(m)=X(m) + x(k)*exp(-j*2*pi*(m-1)*(k-1)/N);
end
end
所以我們只要把程式改成像上面這樣就可以用fft來取代了
像這樣：
X=zeros(N,length(t)); % N=2*Q+1
for n=Q+1:length(t)-Q
p=n-Q:n+Q;
x1=x(p); % x1為從陣列x中取出某段2*Q+1點的新陣列
for m=1:N
for k=1:N
X(m,n)=X(m,n) + x1(k)*exp(-j*2*pi*(m-1)*(k-1)/N); % N=1/(dt*df)
end
end
end
注意length(f)被N取代，N是做FT的範圍，所以用fft的話，df就不能自訂了
然後第二與第三迴圈的部份就跟fft的迴圈型式一樣，所以用fft來取代：
X=zeros(N,length(t)); % N=2*Q+1
for n=Q+1:length(t)-Q
p=n-Q:n+Q;
x1=x(p); % x1為從陣列x中取出某段2*Q+1點的新陣列
X(:,n) = fft(x1);
end
也可以直接把x(p)丟給fft：
X=zeros(N,length(t)); % N=2*Q+1
for n=Q+1:length(t)-Q
p=n-Q:n+Q;
X(:,n) = fft(x(p));
end
這樣就完成了，很簡單吧
不知道為什麼老師上課時會講得這麼複雜 ^^|
至於Gabor就把x(p)乘一個Gaussian function再丟給fft就好了
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 125.230.23.4