[心得] fftconv()

[lihgong (人生，是一句引用句)]

這是一個 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
推 Skylan:推一下 ^^