Быстрое преобразование Фурье (FFT) и увеличение частоты

[mixer94]

Привет форумчане, небольшой вопрос к знатокам, я не силён в математике и в фурье не в силах разбираться. Однако нужен алгоритм питч-шифтинга. Знаю, в интернете полным полно, но они все требовательны к ресурсам а мне в микроконтроллер запихивать надо будет.

Смещение частот (на октаву) по простой формуле:

for i:=1023 downto 0 do
  begin
     arr[2*i]:=arr[trunc(i/2)*2]; //чётные индексы - действительные 
     arr[2*i+1]:=arr[trunc(i/2)*2+1]; //нечётные индексы - мнимые
  end;

Результаты своих экспериментов приложил, обратное преобразование работает прекрасно, но при смещении частот к центру уходит в ноль. Подскажите где накосячил?

[raxp]

...во-первых - к какому центру смещение частот? Куда смотреть? Оси к тому же не подписаны и без значений.
...во-вторых - "различия с исходным красным", где красное? Есть черное и зеленое, красного не видать, куда смотреть?
...в-третьих - какой длительности сэмпл анализируется? Т.е. какой кусок подвергается БПФ, сдвигу гармоник и обратному БПФ?

При использовании смещения частотного спектра требуется проведение дополнительной фильтрации.

[mixer94]

1. простое умножение на 2 выходит. знаю, алгоритм не верный, вот интересуюсь как его сделать чтобы реально менял тональность.
2. Оптимизация jpeg взяла своё...


3. 1024 выборки

[raxp]

Вот, к примеру у меня: генерация исходного тона 500 Гц, производится прямой БПФ (красное), потом сдвигаем исходные гармоники в два раза вниз (2*i) - передаем в обратный БПФ и делаем прямой БПФ для визуализации (зеленое, основной тон 250 Гц +-), потом сдвигаем исходные гармоники в два раза вверх (i div 2) - передаем в обратный БПФ и делаем прямой БПФ для визуализации (синее, основной тон 1000 Гц +-):


(подставка-смещение спектров по вертикали введены для удобства визуализации)

Смотреть видео: http://youtu.be/HnX9Se9bAjU

  // прямое БПФ
  fftio(fDataBuf, 1024, true);

  // сдвиг частот вниз в два раза
  for i:= 0 to length(fDataBuf) -1 do begin
   setlength(fDataOutBuf1, length(fDataOutBuf1)+1);
   if 2*i< length(fDataBuf) then
    fDataOutBuf1[length(fDataOutBuf1)-1]:= fDataBuf[2*i]
     else fDataOutBuf1[length(fDataOutBuf1)-1]:= 0
  end;

  // сдвиг частот вверх в два раза
  for i:= 2 to length(fDataBuf) -1 do begin
   setlength(fDataOutBuf2, length(fDataOutBuf2)+1);
   fDataOutBuf2[length(fDataOutBuf2)-1]:= fDataBuf[i div 2]
  end;

  // обратное FFT для октавы вниз
  fftio(fDataOutBuf1, 1024, false);
  // обратное FFT для октавы вверх
  fftio(fDataOutBuf2, 1024, false);

  // прямое БПФ
  fftio(fDataOutBuf1, 1024, true);
  fftio(fDataOutBuf2, 1024, true);
...

Оптимизация jpeg взяла своё...

используйте нормальные хостинги картинок и выкладывайте в png.

[mixer94]

Спасибо за подробный ответ, возникла парочка вопросов.
Одну ошибку в программе исправил - недорисовывался спектр, он на самом деле зеркальный. Вот и возник вопрос, что делать со второй половиной перед отправкой на обратное преобразование?

Далее у меня два массива чисел после прямого преобразования - действительные и мнимые если не ошибаюсь, в вашем коде этой пары не увидел, поэтому хотелось бы узнать какие из них в какой обработке нуждаются..(перепробовал все варианты однако правильного результата так и не увидел). Ваша библиотека конечно понравилась, но применить её не будет возможности.

[raxp]

...зеркалку отсекать.
...работа с квадратурами у меня для удобства скрыта от глаз.

Ваша библиотека конечно понравилась, но применить её не будет возможности.

в контексте с МК это понятно , главное сам принцип шифтинга.

[mixer94]

не могли бы вы привести эти самые формулы квадратур при прямом и обратном?

[raxp]

...в блоге программистов есть статья.

[mixer94]

Нашёл...всёравно не работает....выкладываю последний результат и весь код, прошу помочь.

procedure FFT(var a : Array of single;
     nn : Integer;
     InverseFFT : Boolean);
var
    ii, jj, n, mmax, m, j, istep, i, isign : Integer;
    wtemp, wr, wpr, wpi, wi, theta, tempr, tempi : Double;
begin
    if InverseFFT then isign := -1
    else isign := 1;
    n := 2*nn; j := 1; ii:=1;
    while ii <= nn do
    begin
        i := 2*ii-1;
        if j>i then
        begin
            tempr := a[j-1];
            tempi := a[j];
            a[j-1] := a[i-1];
            a[j] := a[i];
            a[i-1] := tempr;
            a[i] := tempi;
        end;
        m := n div 2;
        while (m>=2) and (j>m) do
        begin
            j := j-m;
            m := m div 2;
        end;
        j := j+m;
        Inc(ii);
    end;
    mmax := 2;
    while n>mmax do
    begin
        istep := 2*mmax;
        theta := 2*Pi/(isign*mmax);
        wpr := -2.0*sqr(sin(0.5*theta));
        wpi := sin(theta);
        wr := 1.0;
        wi := 0.0;
        ii:=1;
        while ii<=mmax div 2 do
        begin
            m := 2*ii-1;
            jj:=0;
            while jj<=(n-m) div istep do
            begin
                i := m+jj*istep;
                j := i+mmax;
                tempr := wr*a[j-1]-wi*a[j];
                tempi := wr*a[j]+wi*a[j-1];
                a[j-1] := a[i-1]-tempr;
                a[j] := a[i]-tempi;
                a[i-1] := a[i-1]+tempr;
                a[i] := a[i]+tempi;
                Inc(jj);
            end;
            wtemp := wr;
            wr := wr*wpr-wi*wpi+wr;
            wi := wi*wpr+wtemp*wpi+wi;
            Inc(ii);
        end;
        mmax := istep;
    end;
    if InverseFFT then
    begin
        I:=1;
        while I<=2*nn do
        begin
            a[I-1] := a[I-1]/nn;
            Inc(I);
        end;
    end;
end;


procedure draw();
   var i: integer;   
   var k: integer;
begin
  with form1 do begin

  q:=1000;

  for i:=0 to 1023 do
  begin
     arr[2*i]:=sin((q)*i/11024);
     arr[2*i+1]:=0;
     image1.Canvas.LineTo(i,trunc(20*arr[2*i])+50);
  end;

  //прямое преобразование
  fft(arr,1024,false);

  for i:=0 to 1023 do
  begin
     image3.Canvas.MoveTo(i,320);
     image3.Canvas.LineTo(i,230-trunc(sqrt(sqr(arr[2*i])+sqr(arr[2*i+1]))));
  end;

  for i:=0 to 1023 do
  begin
    //увеличение на октаву
    //такая сложность для предотвращения перемешивания Re и Im
    //получение чётного числа + остаток от деления в случае нечётности
      arr2[i]:=arr[(i div 2) div 2 * 2 + (i mod 2)];
      arr2[i+1024]:=0;
  end;

  for i:=0 to 1023 do
  begin
     image5.Canvas.MoveTo(i,330);
     image5.Canvas.LineTo(i,230-trunc(sqrt(sqr(arr2[2*i])+sqr(arr2[2*i+1]))));
  end;

  //обратное для массива с исходными данными
  fft(arr,1024,true);

  //обратное для массива с изменёнными данными
  fft(arr2,1024,true);

  for i:=0 to 1023 do
  begin
     form1.image2.Canvas.LineTo(i,trunc(20*arr[2*i])+50);
  end;

  for i:=0 to 1023 do
  begin
     image4.Canvas.LineTo(i,trunc(20*arr2[2*i])+50);
  end;

  for i:=0 to 1023 do begin
    for k:=0 to 250 do begin
      if Image1.Canvas.Pixels[i,k] <> Image2.Canvas.Pixels[i,k] then
        Image2.Canvas.Pixels[i,k] := clRed;
    end;
  end;
  end;

end;

[raxp]

..на основании каких соображений вы заполняете отсчетами синуса перед БПФ только каждым вторым?

for i:=0 to 1023 do
  begin
     arr[2*i]:=sin((q)*i/11024);
     arr[2*i+1]:=0;
     image1.Canvas.LineTo(i,trunc(20*arr[2*i])+50);
  end;

  //прямое преобразование
  fft(arr,1024,false);

Это вы не со спектром имеете дело, перетасовывать надо после БПФ. И опять же, в последнем случае на выходе вы уже не работаете с квадратурами:

for i:=0 to 1023 do
  begin
    //увеличение на октаву
    //такая сложность для предотвращения перемешивания Re и Im
    //получение чётного числа + остаток от деления в случае нечётности
      arr2[i]:=arr[(i div 2) div 2 * 2 + (i mod 2)];
      arr2[i+1024]:=0;
  end;

тут их нет. См. мой код сдвига выборки выше.