[心得] fftconv()

看板MATLAB作者lihgong (人生，是一句引用句)時間17年前 (2007/07/07 08:58)推噓2(2推 0噓 0→)

留言2則, 2人參與討論串1/1

這是一個 Matlab 的 conv() 的替代品我的程式裡, 需要算大量的 convolution 因為 Matlab 內建的 conv(), 碰到比較長的信號, 速度整個就慢下來為了加速這部分的運算, 我用 DSP 的 overlap-add 的作法, 寫了 fftconv() ---- 如果輸入信號 x 和 response h 都是 vector 那下面這個 function 的表現和 Matlab 完全一樣做同樣的運算, conv() 和 fftconv() 的差異僅止於 machine precision 在做 overlap-add 的時候, 必須把信號切段處理每一段的長如果是 2^n, 那麼 FFT 可以算得飛快所以函數會自動根據輸入長度, 自動選擇適當長度的 FFT-window 切段的次數越少, 速度越快, 不過這需要比較大的記憶體程式會選擇最大可接受的 FFT-window 函數裡, 標藍色的地方要簡單的修改一下, 2^21表示 FFT-window 最大 2 Mega 2^21 是在配備 1G 記憶體的電腦上 tune 出來的設小一點不會出事, 頂多算慢一點點, 設太大用光記憶體, 那就一點都沒得加速了 XD ---- function y = fftconv(x, h) %% 先把輸入整型一下 (to column vector) [dim_x, dim_y] = size(x); if(dim_y > dim_x) x = x'; end [dim_x, dim_y] = size(h); if(dim_y > dim_x) h = h'; end if(length(h) > length(x)) tmp = h; h = x; x = tmp; clear tmp; end % parameter of conv-FFT len_y = length(h) + length(x) - 1; y = zeros(len_y, 1); fft_size = 2^ceil(log2(len_y)); if(fft_size > 2^21) fft_size = 2^21; if(fft_size < length(h)) fprintf('Too Large input...\n'); end end processing_x_block_size = fft_size +1 - length(h); h_freq = fft(h, fft_size); clear h; %% 算 FFT 喔~~~ x_index_begin = 1; while x_index_begin < length(x) x_index_end = x_index_begin + processing_x_block_size - 1; if(x_index_end > length(x)) x_index_end = length(x); end tmp = ifft( ... fft(x(x_index_begin:x_index_end), fft_size) .* ... h_freq); % 把 output 加上去 y_index_begin = x_index_begin; y_index_end = x_index_begin+length(tmp)-1; if(y_index_end > length(y)) y_index_end = length(y); end len_valid_output = y_index_end - y_index_begin + 1; y(y_index_begin:y_index_end) = ... y(y_index_begin:y_index_end) + tmp(1:len_valid_output); x_index_begin = x_index_begin + processing_x_block_size; end -- ※ 發信站: 批踢踢實業坊(ptt.cc) ◆ From: 140.113.236.184