参考数字信号处理教程第四版程佩青著第四章FFT这里直接给出matlab函数性能不保证最优注意此函数只能处理混合基fft即输入信号x的长度不能是素数不能是2次幂整数function X mixedRadixFFT(x) % multiBasesFFT3 - 任意混合基FFT运算 % 输入: x - 任意长度的数据行向量 % 输出: X - 频域结果列向量 validateattributes(x, {double},{row, nonempty}); N numel(x); % 输入信号点数 factorSeq factor(N); % 找出所有的基 r getr1r2torLStyle(factorSeq); % 整合所有为2的基输出符合要求的r0r1...rL-1形式 L numel(r); % 基的数量 % 获取x(n)的dft索引顺序和最后输出的X(k)的倒位序 [nCoef, kCoef] getCoef(r); % 计算n和k的系数 n 0; k 0; % 获取nL-1,...,n0用于计算旋转因子 nRepsPerMatForEach cumprod(r); nArr zeros(N, L); for i 1 : L % 计算nArr nL-1,...,n0 ri r(end - (i - 1)); nRepsPerElem cumprod(r(end - (i - 1) : end)); nRepsPerElem N / nRepsPerElem(end); ni repelem((0 : ri - 1), nRepsPerElem); nRepsPerMat N / nRepsPerMatForEach(end - (i - 1)); nArr(:, i) repmat(ni, [nRepsPerMat 1]); % 获取x(n)的正序n r0r1...rL-2*nL-1 r0r1...rL-3*nL-2 ... r0r1*n2 r0*n1 n0 n n(:) nCoef(i) * (0 : r(end - (i-1)) - 1); % 获取X(k)的倒位序k r1r2...rL-1*n0 r2r3...rL-1*n1 ... rL-2rL-1*nL-3 rL-1*rL-2 nL-1 k k(:) kCoef(i) * (0 : r(i) - 1); end r r(end : -1 : 1); % 反转r便于索引 ind reshape(n(:), r); % 每一个维度对应每个基的dft索引 indCell cell(L, 1); d (1 : L); % matlab角度的维度索引用于permute % 混合基FFT计算 Xtmp x.; for iD 1 : L % 索引每个基 % 获取当前基的dft索引转换为2维矩阵方便观察和计算 dtmp d(d ~ iD); dftInd reshape(permute(ind, [iD dtmp]), r(iD), []); indCell{iD} dftInd; % 记录 % 对每一个维度进行dft W_N exp(-1i * 2 * pi / r(iD) * (0 : r(iD) - 1) .* (0 : r(iD) - 1)); Xtmp(dftInd 1) Xtmp(dftInd 1) * W_N; % 如果当前并非最后一个基则需要点乘上旋转因子 if iD ~ L % nRepsPerMatForEach(end - (iD - 1))为当前旋转因子的N twiddleFactor exp(-1i * 2 * pi / nRepsPerMatForEach(end - (iD - 1)) * sum(nCoef(iD 1 : end) .* nArr(:, iD 1 : end), 2) .* nArr(:, iD)); Xtmp Xtmp .* twiddleFactor; end end % 整序输出 X zeros(N, 1); X(k(:) 1) Xtmp; end %% Local functions function [nCoef, kCoef] getCoef(r1r2torL) % getCoef - 获取n和k的系数也就是(n)10和(bar{n})10的系数 % 输入r1r2torL - 基行向量或列向量 % 输出nCoef - n的系数 % 输出kCoef - k的系数 n numel(r1r2torL); res cumprod(r1r2torL); res(end) []; nCoef [res(end : -1 : 1) 1]; kCoef ones(1, n); for i 2 : n res cumprod(r1r2torL(i : end)); kCoef(i - 1) res(end); end end function r1r2torL getr1r2torLStyle(factorSeq) % getr1r2torLStyle - 输出r1r2...rL形式的混合基其中r1整合了所有基2 % 输入factorSeq - 将N分解出所有的素因子 % 输出r1r2torL - r1r2...rL形式的混合基 prodRes cumprod(factorSeq(factorSeq 2)); if isempty(prodRes) r1 []; else r1 prodRes(end); end r2torL factorSeq(factorSeq ~ 2); r1r2torL [r1 r2torL]; end
