我修复了这个错误。此处的工作代码
Rader DFT算法使用GNU Octave实现(例如，长度为11(。我用过这篇维基百科文章。获得的值是正确的，但它们被错误地重新索引。我不明白错误在哪里？
upd。的Add函数可查找组中最小的生成器。
Fin = [1,2,3,4,5,6,7,8,9,10,11];
nfft = columns(Fin);
snfft = nfft - 1;
function [m one_per] = find_gen(p)
m = 1;
a = factor(p);
if (length(a) > 1)
m = p + 1;
one_per = [];
endif
if (m == 1)
finished = 0;
for cur = 2:p-2
not_yet = 0;
test = cur;
single_per = [1 cur];
for k = 2:p-2
test = test * cur;
test = mod(test,p);
single_per = [single_per test];
if (test == 1)
not_yet = 1;
endif
endfor
if (not_yet == 0)
m = cur;
one_per = single_per;
finished = 1;
break;
endif
endfor
endif
endfunction
q = find_gen(nfft)
p = mod((q^(nfft-2)),nfft)
Tq_idx = [];
Tq = [];
for k = 0 : snfft-1
A = mod(q^k, nfft);
Tq_idx = [A Tq_idx];
Tq = [Fin(A+1) Tq];
endfor
Tq_idx, Tq
Tp_idx = [];
Tp = [];
for k = 0 : snfft-1
A = mod(p^k, nfft);
Tp_idx = [A Tp_idx];
Tp = [Fin(A+1) Tp];
endfor
Tp_idx, Tp
Twp = [];
for k = 1 : snfft
ecpx = complex(cos(-2*pi*Tp_idx(k) / nfft),sin(-2*pi*Tp_idx(k) / nfft));
Twp = [Twp ecpx];
endfor
Tq_fft = fft(Tq);
Twp_fft = fft(Twp);
Tm_fft = Tq_fft .* Twp_fft;
Tm_ffti = ifft(Tm_fft);
Tm_ffti
Res = [Fin(1)];
for k = 1 : snfft
Res(1) += Fin(k+1);
Res = [Res (Tm_ffti(Tp_idx(k)) + Fin(1))];
endfor
Res
Fbest = fft(Fin)
Fdiff = Fbest .- Res
ResI = ifft(Res)
结果
Res =
Columns 1 through 3:
66.0000 + 0.0000i -5.5000 - 4.7658i -5.5000 - 18.7313i
Columns 4 through 6:
-5.5000 - 0.7908i -5.5000 - 8.5582i -5.5000 + 8.5582i
Columns 7 through 9:
-5.5000 + 18.7313i -5.5000 + 4.7658i -5.5000 + 0.7908i
Columns 10 and 11:
-5.5000 - 2.5118i -5.5000 + 2.5118i
使用GNU Octave fft((内部函数作为标准
Fbest =
Columns 1 through 3:
66.0000 + 0.0000i -5.5000 + 18.7313i -5.5000 + 8.5582i
Columns 4 through 6:
-5.5000 + 4.7658i -5.5000 + 2.5118i -5.5000 + 0.7908i
Columns 7 through 9:
-5.5000 - 0.7908i -5.5000 - 2.5118i -5.5000 - 4.7658i
Columns 10 and 11:
-5.5000 - 8.5582i -5.5000 - 18.7313i
function [m one_per] = find_gen(p)
m = 1;
a = factor(p);
if (length(a) > 1)
m = p + 1;
one_per = [];
endif
if (m == 1)
finished = 0;
for cur = 2:p-2
not_yet = 0;
test = cur;
single_per = [1 cur];
for k = 2:p-2
test = test * cur;
test = mod(test,p);
single_per = [single_per test];
if (test == 1)
not_yet = 1;
endif
endfor
if (not_yet == 0)
m = cur;
one_per = single_per;
finished = 1;
break;
endif
endfor
endif
endfunction
function [Fout] = fast_conv (F1,F2, seq_lenght)
F1_fft = fft(F1);
F2_fft = fft(F2);
Fm_fft = F1_fft .* F2_fft;
Fout = ifft(Fm_fft);
endfunction
function [Res] = rader_algo (Fin)
nfft = columns(Fin);
snfft = nfft - 1;
q = find_gen(nfft)
p = mod((q^(nfft-2)),nfft)
Tq_idx = []; Tq = [];
for k = 0 : snfft-1
A = mod(q^k, nfft);
Tq_idx = [Tq_idx A];
Tq = [Tq Fin(A+1)];
endfor
Tq_idx, Tq
Tp_idx = []; Tp = [];
for k = 0 : snfft-1
A = mod(p^k, nfft);
Tp_idx = [Tp_idx A];
Tp = [Tp Fin(A+1)];
endfor
Tp_idx, Tp
Twp = [];
for k = 1 : snfft
ecpx = complex(cos(-2*pi*Tp_idx(k) / nfft),sin(-2*pi*Tp_idx(k) / nfft));
Twp = [Twp ecpx];
endfor
Twp
Tm_ffti = fast_conv(Tq, Twp, nfft);
Res = zeros(1, nfft);
Res(1) = Fin(1);
for k = 1 : snfft
Res(1) += Fin(k+1);
Res(Tp_idx(k)+1) = Tm_ffti(k) + Fin(1);
endfor
endfunction
Fin = [1,2,3,4,5];
Res = rader_algo (Fin)
% === VERIFY ===
Fbest = fft(Fin)
Fdiff = Fbest .- Res
ResI = ifft(Res)