Частотный преобразователь с использованием FFT

Question

Я работаю на частотном переключателе с использованием примитивного алгоритма FFT, поставляемого Rosetta Code. Я понимаю, что для сдвига частоты сигнала образцов, один применяет БПФ к исходному звуку, умножает частоту каждой результирующей синусоидальной частоты на частотный сдвиг (определяемый пользователем), а затем снова объединяет синусоидальные волны. Когда я запускаю свой алгоритм, выход имеет крайне низкое качество, как будто в алгоритме не было достаточного количества синусоидальных волн, чтобы вначале правильно воспроизвести сигнал. Алгоритм реализован в классе в файле заголовка и называется (правильно) в другом месте.

#include <complex> 
#include <valarray> 

typedef std::complex<double> Complex; 
typedef std::valarray<Complex> CArray; 

class FrequencyShifter { 
float sampleRate; 
public: 
    FrequencyShifter() { 

    } 
    void setSampleRate(float inSampleRate) { 
     sampleRate = inSampleRate; 
    } 
    double abs(double in0) { 
     if (in0>=0) return in0; 
     else return -in0; 
    } 
    void fft(CArray& x) 
    { 
     const size_t N = x.size(); 
     if (N <= 1) return; 

     // divide 
     CArray even = x[std::slice(0, N/2, 2)]; 
     CArray odd = x[std::slice(1, N/2, 2)]; 

     // conquer 
     fft(even); 
     fft(odd); 

     // combine 
     for (size_t k = 0; k < N/2; ++k) 
     { 
      Complex t = std::polar(1.0, -2 * PI * k/N) * odd[k]; 
      x[k ] = even[k] + t; 
      x[k+N/2] = even[k] - t; 
     } 
    } 
    double convertToReal(double im, double re) { 
     return sqrt(abs(im*im - re*re)); 
    } 
    void processBlock(float *inBlock, const int inFramesToProcess, float scale) { 
     //inFramesToProcess is the amount of samples in inBlock 
     Complex *copy = new Complex[inFramesToProcess]; 
     for (int frame = 0; frame<inFramesToProcess; frame++) { 
      copy[frame] = Complex((double)inBlock[frame], 0.0); 
     } 
     CArray data(copy, inFramesToProcess); 
     fft(data); 
     const float freqoffsets = sampleRate/inFramesToProcess; 
     for (float x = 0; x<data.size()/2; x++) { 
      for (float frame = 0; frame<inFramesToProcess; frame++) { 
       inBlock[(int)frame] = (float)(convertToReal(data[(int)x].imag(), data[(int)x].real())*sin(freqoffsets*x*frame*scale)); 
      } 
     } 
    } 
}; 

Я предполагаю, что часть проблемы заключается в том, что я только в том числе sampleRate/inFramesToProcess частоты для синусоиды, чтобы покрыть. Посылал бы большие аудиофайлы (таким образом, больше *inBlock s и inFramesToProcess s) сделать звук менее зернистым? Как бы это сделать без, просто изменяя значения или длины аргументов?

Answer 1

Вот обновленная версия processBlock с количеством настроек, необходимых для реализации частотного сдвига, который я опишу ниже:

void processBlock(float *inBlock, const int inFramesToProcess, float scale) { 
    //inFramesToProcess is the amount of samples in inBlock 
    Complex *copy = new Complex[inFramesToProcess]; 
    for (int frame = 0; frame<inFramesToProcess; frame++) { 
     copy[frame] = Complex((double)inBlock[frame], 0.0); 
    } 
    CArray data(copy, inFramesToProcess); 
    fft(data); 
    const float freqoffsets = 2.0*PI/inFramesToProcess; 
    const float normfactor = 2.0/inFramesToProcess; 
    for (int frame = 0; frame<inFramesToProcess; frame++) { 
     inBlock[frame] = 0.5*data[0].real(); 
     for (int x = 1; x<data.size()/2; x++) { 
      float arg = freqoffsets*x*frame*scale; 
      inBlock[frame] += data[x].real()*cos(arg) - data[x].imag()*sin(arg); 
     } 
     inBlock[frame] *= normfactor; 
    } 
} 

Вывод

Спектр Вы получаете от быстрого преобразования Фурье является комплекснозначным, что можно рассматривать как обеспечение представления ваших сигналов в терминах синусоидальных и косинусных волн. Восстанавливая временной области сигнала может быть сделано с помощью обратного преобразования, который будет задается соотношением:

Воспользовавшись симметрией частотного спектра, то это может быть выражено как:

или что то же самое:

Как вы могли заметить, этот термин по индексу 0 и N/2 - частные случаи с чисто вещественными коэффициентами в частотной области. Для простоты, предполагая, что спектр не проходит до N/2, вы можете отказаться от этого термина N/2 и по-прежнему получать разумное приближение. Для других условий вы получите вклад, который может быть реализован в виде

normfactor = 2.0/inFramesToProcess; 
normfactor*(data[x].real()*cos(arg) - data[x].imag()*sin(arg)) 

Вы бы, конечно, нужно добавить все эти взносы в буфер окончательной inBlock[frame], а не просто перезапись предыдущих результатов:

inBlock[frame] += normfactor*(data[x].real()*cos(arg) - data[x].imag()*sin(arg)); 
//    ^^ 

Обратите внимание, что нормализация может быть выполнена в конечном результате после цикла, чтобы уменьшить количество умножений.При этом, мы должны обратить особое внимание на термин DC с индексом 0 (который имеет коэффициент 1/N вместо 2/N):

inBlock[frame] = 0.5*data[0].real(); 
for (int x = 1; x<data.size()/2; x++) { 
    float arg = freqoffsets*x*frame*scale; 
    inBlock[frame] += data[x].real()*cos(arg) - data[x].imag()*sin(arg); 
} 
inBlock[frame] *= normfactor; 

Наконец, при генерации тональных сигналов, фазовый аргумент arg в sin и cos должен иметь вид 2*pi*k*n/inFramesToProcess (до применения фактора scale), где n - индекс выборки во временной области, а k - индекс частотной области. Конечным результатом является то, что вычисленный приращение частоты freqoffsets действительно должен быть 2.0*PI/inFramesToProcess.

Примечания

• Алгоритм БПФ работает на предположении, что ваш основной сигнал во временной области является периодическим с периодом вашей длины блока. В результате могут быть слышимые разрывы между блоками.
• Будущим читателям следует знать, что это не меняет спектр на постоянное количество, а скорее хлюпает или расширяет спектр по мере того, как частоты масштабируются мультипликативным фактором. Например, сигнал, который включает в себя компоненты 100-200 Гц, может быть сжат до 75-150 Гц в 0,75 раза. Обратите внимание, что нижний предел был сдвинут на 25 Гц, а верхний предел был понижен на 50 Гц.