function x = fft_2n(x)
% 計算長度為2n的實序列的fft
% john g.proakis,數字訊號處理,第四版
% define variables:
% x --- input real vector of length 2*n(power of 2)
% x --- output vector,fft result
% record of revisions:
% date programmer description of change
% 2014/0421 gou dh original code
% test
clear,clc
x = randn(64,1);
x = x(:);
nx = length(x);
t = reshape(x,2,nx/2);
x1 = t(1,:)';
x2 = t(2,:)'; % 時域抽取,奇偶分離
[x1,x2] = fft_2_n(x1,x2,nx/2);
k = [0:nx/2-1]';
w2nk = exp(-j*2*pi*k/(nx));
t = x2.*w2nk;
g1 = x1 + t;
g2 = x1 - t;
x = [g1;g2];
% test
sum(abs(x - fft(x,nx)))
return;
% subfunction
function [x1,x2] = fft_2_n(x1,x2,n)
% 計算2個n點實序列的fft
if ~isreal(x1) | ~isreal(x2)
return;
endx = x1 + j*x2;
x = fft_dif2(x,-1);
t = [x(1); flipud(x(2:end))];
t = conj(t);
x1 = (x + t)/2;
x2 = (x - t)/(2*j);
% subfunction
function [x]=fft_dif2(x,f)
% 頻域抽取基2-fft(基2-dif-fft)
nx = length(x);
m = fix(log2(nx));
if 2^m < nx, m = m + 1; end % 級數m
x = [x(:);zeros(2^m-nx,1)];
n = length(x);
c = j*pi*f;
for m = 0:1:m-1 % 級迴圈序號
b = 2^(m-1-m); % 蝶形因子序號間距=蝶形因子個數=2^(m-1-m)
for r = 0:1:b-1 % 蝶形因子迴圈序號
wnp = exp(c*r/b);
for k = r:2*b:n-1
tmp1 = x(k+1);
tmp2 = x(k+b+1);
x(k+1) = tmp1+tmp2;
x(k+b+1) = (tmp1-tmp2)*wnp;
endend
end% 倒位序
% x = x;
lh = n/2;
j = lh;
n1 = n-2;
for i = 1:1:n1
if i < j
t = x(i+1);
x(i+1) = x(j+1);
x(j+1) = t;
endk = lh;
while j >= k
j = j-k;
k = k/2;
endj = j+k;
endif f == 1
x = x/n;
endx = x;
收藏於 2014-04-21
DIT基2點fft的程式設計實現
找了很多資料，發現這個最有用，思路簡單清晰，結合自己的理解貼出來。徹底弄懂基2點流程之後，基4也就不難。下一次貼出基4fft程式設計實現。蝶形公式 x k x k x k b w pn x k b x k x k b w pn 其中w pn cos 2 p n jsin 2 p n 設 x k b ...
找大於等於乙個數的最小的2 n
最近看hashmap原始碼時，發現給定初始capacity計算threshold的過程很巧妙。1 static final int tablesizefor int cap 這裡是實現了找大於等於cap的最小2 n。cap為int型別，長度32位。對於乙個正數，找大於該數的最小的2 n,都可以採用這...
求乙個數字大於並最接近的2 N
無符號右移 無符號右移運算子 同右移，但是結果全變正數。或 或運算 二進位制中只要乙個為1就為1 在hashmap原始碼中有相關操作，直接分析原始碼如下 返回給定目標容量的2倍冪。將我們傳入的容量設定為大於並最接近的2 n 補位，將原本為0的空位填補為1，最後加1時，最高有效位進1，其餘變為0,如此...