这是一个 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