Неправильная индексация в алгоритме Rader (реализация GNU Octave)

Question

Алгоритм Rader DFT реализован с использованием GNU Octave (например, длина 11). Я использовал эту статью в Википедии. Полученные значения верны, но они неправильно переиндексированы. Не могу понять, где ошибка?

upd. Функция добавления для находит самый маленький генератор группы.

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

Answer 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

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)