看板 DJJ_Course
作者 標題 [TFW]作業三的Gabor Transform使用CZT
時間 2007年12月03日 Mon. AM 01:34:05
我終於想起來了
在第二題要使用CZT的Garbor Tranform來做ringin.wav為什麼會出問題
原來問題在於使用CZT時，有個限制是df要與dt相同
也就是頻率解析度與時間解析度要一樣
至於為何？我想是因為CZT的運算中沒有辦法自訂df的大小吧
如果只是做一般的訊號也沒什麼關係
反正這兩個解析度就愈高愈好嘛
不過用在ringin.wav就不行了
因為這個訊號的取樣頻率太高了，fs=11025
所以dt=1/fs≈0.00009
而df又必需等於dt
所以頻率軸必需取 f=-fs/2:df:fs/2
頻率軸的長度M會高達 11025*11025
太大了，會造成out of memory
於是我想到了一個方法，就是自己把取樣頻率降低
例如 fs=fs/100=110.25
dt=1/fs=0.0091
頻率軸的長度M即降為 110.25*110.25≈10156
而對訊號的影響就是頻寬縮小100倍，時間長度變為原本的100倍
只要在顯示的時候修改一下t軸與f軸即可
而且這樣順便解決了一個問題
Gabor Transform中用到的Gaussian function的寬度約1秒
但訊號ringin.wav的長度根本不到1秒
所以求得的時頻分佈圖會看不出頻率對時間的變化
但現在時間長度變成原本的100倍後，就沒有這個問題了
只是這樣做，還是會出現out of memory
所以ringin.wav這個訊號我只取9千多點中的前面的1000點
這樣就可以跑出來了...
程式碼給大家參考一下：
[x,fs]=wavread('ringin.wav');
N = 1000; % 時間軸的長度
x = x(1:N).'; % 訊號x只取前面N點，並轉為row向量
fs=fs/100; % 取樣頻率太高，先降為1/100
dt=1/fs; % 時間軸的解析度、取樣時間
t=(0:N-1)*dt; % 時間軸
% 產生頻率軸f
df=dt; % 頻率軸的解析度必需和時間軸的解析度相同
f=-fs/2:df:fs/2; % 頻率軸取原點附近一個fs的區間
M=length(f); % 頻率軸的長度
t1=ones(N,1)*t; % 代表時間軸t的NxN矩陣
t2=t1.'; % 代表tau的NxN矩陣
gs=exp(-pi*(t2-t1).^2); % Gaussian function
% 產生convolution的輸入訊號x2
x1=(x.'*ones(1,N)).*gs;
x2=x1.*exp(-j*pi*t2.^2);
% 產生convolution的脈衝響應x3, t3為其時間範圍
t3=min(f)-max(t)+(0:M+N-2)*dt;
x3=exp(j*pi*t3.^2).';
% 將x2的column與x3做convolution後，存進x4的column
x4=zeros(M+N*2-2,N);
for n=1:N
x4(:,n)=conv(x2(:,n),x3)*dt;
end
f1=f.'*ones(1,N); % 代表頻率軸的f的NxN矩陣
x4=x4(N:N+M-1,:); % convolution的結果只要中間的M個值
X=x4.*exp(-j*pi*f1.^2);
% 顯示X的時頻分佈圖，將最大值調整至64(白色)
% 時間軸與頻率軸調整回正確的值
image(t/100,f*100,abs(X)*64/max(max(abs(X))));
colormap('Gray'); % 用0(黑)~64(白)的灰階值顯示
執行結果：N = 1000; % 時間軸的長度
x = x(1:N).'; % 訊號x只取前面N點，並轉為row向量
fs=fs/100; % 取樣頻率太高，先降為1/100
dt=1/fs; % 時間軸的解析度、取樣時間
t=(0:N-1)*dt; % 時間軸
% 產生頻率軸f
df=dt; % 頻率軸的解析度必需和時間軸的解析度相同
f=-fs/2:df:fs/2; % 頻率軸取原點附近一個fs的區間
M=length(f); % 頻率軸的長度
t1=ones(N,1)*t; % 代表時間軸t的NxN矩陣
t2=t1.'; % 代表tau的NxN矩陣
gs=exp(-pi*(t2-t1).^2); % Gaussian function
% 產生convolution的輸入訊號x2
x1=(x.'*ones(1,N)).*gs;
x2=x1.*exp(-j*pi*t2.^2);
% 產生convolution的脈衝響應x3, t3為其時間範圍
t3=min(f)-max(t)+(0:M+N-2)*dt;
x3=exp(j*pi*t3.^2).';
% 將x2的column與x3做convolution後，存進x4的column
x4=zeros(M+N*2-2,N);
for n=1:N
x4(:,n)=conv(x2(:,n),x3)*dt;
end
f1=f.'*ones(1,N); % 代表頻率軸的f的NxN矩陣
x4=x4(N:N+M-1,:); % convolution的結果只要中間的M個值
X=x4.*exp(-j*pi*f1.^2);
% 顯示X的時頻分佈圖，將最大值調整至64(白色)
% 時間軸與頻率軸調整回正確的值
image(t/100,f*100,abs(X)*64/max(max(abs(X))));
colormap('Gray'); % 用0(黑)~64(白)的灰階值顯示
![[圖]](http://knuckles.twbbs.org/photo/temp/ringin.jpg)
--
※ 來源: 台大電信 DISP 實驗室 (http://disp.twbbs.org)
※ 作者: Knuckles 來自: 140.112.175.130 時間: 2007-12-03 01:34:06
※ 編輯: Knuckles 來自: 140.112.175.130 時間: 2007-12-03 01:37:35
推 dog: 歐歐~~~只要做正頻就行囉...可省一半 >>140.112.175.132 12-03 10:53
推 cnnegro: 神人阿德!! >>59.117.193.60 12-03 23:01
推 nick: 福德正神! >>140.112.175.133 12-04 19:08
推 nick: 如果再加上把Gaunssian function作scaling(變窄)應該更好 >>140.112.175.133 12-06 10:55
※ 編輯: Knuckles 來自: 140.112.175.128 時間: 2008-10-13 04:18:21
※ 編輯: Knuckles 時間: 2010-10-23 05:02:06 來自: 111-248-0-184.dynamic.hinet.net
※ 看板: DJJ_Course 文章推薦值: 4 目前人氣: 0 累積人氣: 844
回列表(←)
分享