Реализация собственного алгоритма быстрого преобразования Фурье не работает - PullRequest
Новая
       6

Реализация собственного алгоритма быстрого преобразования Фурье не работает

0 голосов
/ 06 декабря 2018

Я нашел и настроил алгоритм быстрого преобразования Фурье.Алгоритм приведен ниже: 

function X = myFFT(x)
   %only works if N = 2^k
   N = numel(x);
   xp = x(1:2:end);
   xpp = x(2:2:end);
   if N>=8
      Xp = myFFT(xp);
      Xpp = myFFT(xpp);
      X = zeros(N,1);
      Wn = exp(-1i*2*pi.*((0:N/2-1)')/N);
      tmp = Wn .* Xpp;
      X = [(Xp + tmp);(Xp -tmp)];
   else
     switch N
       case 2
          X = [1 1;1 -1]*x;
       case 4
          X = [1 0 1 0; 0 1 0 -1i; 1 0 -1 0;0 1 0 1i]*[1 0 1 0;1 0 -1 0;0 1 0 1;0 1 0 -1]*x;
       otherwise
          error('N not correct.');
     end
   end
end

Приведенный выше алгоритм хранится в файле .m с именем myFFT.

Я хочу использовать код для двух последующих сигналов: h1 = sin (t2)и для времени t2 = 0: 1: 255;

Я вспоминаю функцию myFFT в другом файле .m следующим образом: 

x=h1;
X11 = myFFT(x);

Но это выдает мне эту ошибку: 

Error using  * 
Inner matrix dimensions must agree.

Error in myFFT (line 18)
           X = [1 0 1 0; 0 1 0 -1i; 1 0 -1 0;0 1 0 1i]*[1 0 1 0;1 0 -1 0;0 1 0 1;0 1 0 -1]*x;

Error in myFFT (line 7)
    Xp = myFFT(xp);

Error in myFFT (line 7)
    Xp = myFFT(xp);

Error in myFFT (line 7)
    Xp = myFFT(xp);

Error in myFFT (line 7)
    Xp = myFFT(xp);

Error in myFFT (line 7)
    Xp = myFFT(xp);

Error in myFFT (line 7)
    Xp = myFFT(xp);

Честно говоря, я понятия не имею, как это исправить.Не могли бы вы мне помочь?Заранее большое спасибо!

1 Ответ

0 голосов
/ 06 декабря 2018

Проблема заключается в умножении матрицы 4х4 на матрицу 1х4.Переменная x должна быть транспонирована.

   function X = myFFT(x)
       %only works if N = 2^k
       N = numel(x);
       xp = x(1:2:end);
       xpp = x(2:2:end);
       if N>=8
          Xp = myFFT(xp);
          Xpp = myFFT(xpp);
          X = zeros(N,1);
          Wn = exp(-1i*2*pi.*((0:N/2-1)')/N);
          tmp = Wn .* Xpp;
          X = [(Xp + tmp);(Xp -tmp)];
       else
         switch N
           case 2
              X = [1 1;1 -1]*x';
           case 4
              X = [1 0 1 0; 0 1 0 -1i; 1 0 -1 0;0 1 0 1i]*[1 0 1 0;1 0 -1 0;0 1 0 1;0 1 0 -1]*x';
           otherwise
              error('N not correct.');
         end
       end
    end

, затем протестируйте его с помощью 

t2 = 0:1:255;
myFFT(sin(t2));
...