Вывод KissFFT из kiss_fftr

Question

Я получаю данные PCM через сокет-соединение в пакетах, содержащих 320 образцов. Частота дискретизации звука составляет 8000 сэмплов в секунду. Я делаю с этим что-то вроде этого:

int size = 160 * 2;//160;
int isinverse = 1;
kiss_fft_scalar zero;
memset(&zero,0,sizeof(zero));
kiss_fft_cpx fft_in[size];
kiss_fft_cpx fft_out[size];
kiss_fft_cpx fft_reconstructed[size];

kiss_fftr_cfg fft = kiss_fftr_alloc(size*2 ,0 ,0,0);
kiss_fftr_cfg ifft = kiss_fftr_alloc(size*2,isinverse,0,0);

for (int i = 0; i < size; i++) {
    fft_in[i].r = zero;
    fft_in[i].i = zero;
    fft_out[i].r = zero;
    fft_out[i].i = zero;
    fft_reconstructed[i].r = zero;
    fft_reconstructed[i].i = zero;
}

// got my data through socket connection

for (int i = 0; i < size; i++) {
     // samples are type of short
     fft_in[i].r = samples[i];
     fft_in[i].i = zero;
     fft_out[i].r = zero;
     fft_out[i].i = zero;
 }

 kiss_fftr(fft, (kiss_fft_scalar*) fft_in, fft_out);
 kiss_fftri(ifft, fft_out, (kiss_fft_scalar*)fft_reconstructed);

 // lets normalize samples
 for (int i = 0; i < size; i++) {
     short* samples = (short*) bufTmp1;
     samples[i] = rint(fft_reconstructed[i].r/(size*2));
 }

После этого я заполняю буферы OpenAL и играю их. Все работает отлично, но я хотел бы сделать некоторую фильтрацию аудио между kiss_fftr а также kiss_fftri, Начальная точка, как я думаю для этого, - преобразовать звук из временной области в частотную, но я не совсем понимаю, какие данные я получаю от kiss_fftr функция. Какая информация хранится в каждом из этих комплексных чисел, какая ее действительная и мнимая часть может рассказать мне о частоте. И я не знаю, какие частоты покрыты (какой диапазон частот) в fft_out - какие индексы соответствуют каким частотам.

Я новичок в области обработки сигналов и преобразований Фурье.

Любая помощь?

Answer 1

Я постараюсь ответить на ваши вопросы напрямую.

// a) the real and imaginary components of the output need to be combined to calculate the amplitude at each frequency. 

float ar,ai,scaling; 

scaling=1.0/(float)size;

// then for each output [i] from the FFT...

ar = fft_out[i].r; 
ai = fft_out[i].i; 
amplitude[i] = 2.0 * sqrtf( ar*ar + ai*ai ) * scaling ;

// b) which index refers to which frequency? This can be calculated as follows. Only the first half of the FFT results are needed (assuming your 8KHz sampling rate) 

for(i=1;i<(size/2);i++) freq = (float)i / (1/8000) / (float)size ; 

// c)  phase (range +/- PI) for each frequency is calculated like this: 

phase[i] = phase = atan2(fft_out[i].i / fft_out[i].r);

Answer 2

Прежде чем перейти обеими ногами к реализации C, ознакомьтесь с цифровыми фильтрами, особенно с фильтрами FIR.

Вы можете спроектировать FIR-фильтр, используя что-то вроде сигнальной панели инструментов GNU Octave. Посмотрите на команду fir1(самая простая), firls или Rememberz. Кроме того, вы можете разработать FIR-фильтр через веб-страницу. Быстрый веб-поиск "онлайн-дизайна фильтра пихты" обнаружил это (я не использовал его, но, похоже, он использует дизайн equiripple, используемый в команде revz или firpm)

Попробуйте сначала внедрить свой фильтр с прямым свертыванием (без БПФ) и посмотрите, приемлема ли скорость - это более простой путь. Если вам нужен подход, основанный на FFT, есть пример реализации overlap-save в файле kissfft/tools/kiss_fastfir.c.

Answer 3

Возможно, вы захотите исследовать быстрое свертывание FFT с использованием алгоритмов добавления с наложением или с перекрытием. Вам нужно будет увеличить длину каждого БПФ на длину импульса желаемого фильтра. Это связано с тем, что (1) свертка FFT/IFFT является круговой, и (2) каждый индекс в результате массива FFT соответствует почти всем частотам (отклик в форме Sinc), а не только одной (даже если в основном около одной), поэтому любой Модификация бина будет просачиваться по всей частотной характеристике (кроме определенных точных периодических частот).