Добрый день Сергей.
Нужна Ваша консультация по следующему вопросу:
Например: у меня на входе гармоника, я выполняю БПФ, затем делаю обратное преобразование Фурье, получаю ту же гармонику на выходе. Мне необходимо получить исходную амплитуду сигнала. При БПФ на каждом шаге бабочки делается масштабирование входных данных для предотвращения переполнения разрядной сетки(деление на 2-т.к. у меня данные с фиксированной точкой 32 разряда, пишу для BlackFin BF518), при 512 точках БПФ у меня получается 9 сдвигов вправо, т.е. деление на 512. При обратном преобразовании опять происходит такое же масштабирование.
Где нужно делать обратное масштабирование?
После прямого БПФ и обратного БПФ или только после обратного БПФ?
Когда я делаю масштабирование после обратного БПФ, то получаю исходную амплитуду сигнала. А если и после прямого и после обратного, то на выходе непонятно что.
вот пример кода
in[NUM_INSAMPLES] - входной сигнал во временной области
out[NUM_INSAMPLES] - выходной массив БПФ(комплексный) в частотной области
out_1[NUM_INSAMPLES] - выходной массив Обратного БПФ(комплексный) в частотной области
block_exponent - тут функции дают число сдвигов = log2(NUM_FFT)=9
NUM_FFT=512
Прямое БПФ: rfft_fr32(in, out, twiddle_table, 1, NUM_FFT, &block_exponent, 1);
Масштабирование for (i=0;i<NUM_INSAMPLES;i++) {
out.re = out_1.re << block_exponent;
out.im = out_1.im << block_exponent;
}
Обратное БПФ: ifft_fr32(out, out_1, twiddle_table, 1, NUM_FFT, &block_exponent, 1);
Масштабирование for (i=0;i<NUM_INSAMPLES;i++) {
out_1.re = out_1.re << block_exponent;
out_1.im = out_1.im << block_exponent;
}
С уважением Михаил,
Ростовская область, г.Таганрог.
БПФ
Re: БПФ
ifft
N-point radix-2 inverse FFT
Synopsis
#include <filter.h>
void ifft_fr16(const complex_fract16 input[],
complex_fract16 output[],
const complex_fract16 twiddle_table[],
int twiddle_stride,
int fft_size,
int *block_exponent,
int scale_method);
void ifft_fr32(const complex_fract32 input[],
complex_fract32 output[],
const complex_fract32 twiddle_table[],
int twiddle_size,
int fft_size,
int *block_exponent,
int scale_method);
Description
The ifft functions transform the frequency domain complex input signal sequence to the time domain by using the radix-2 Fast Fourier Transform.
The size of the input array input and the output array output is fft_size, where fft_size represents the number of points in the FFT. By allocating these arrays in different memory banks, any potential data bank collisions are avoided, thus improving run-time performance. If the input data can be overwritten, the optimum memory usage can be achieved by also specifying the input array as the output array.
The twiddle table is passed in the argument twiddle_table, which must contain at least fft_size/2 twiddle factors. The table is composed of +cosine and -sine coefficients and may be initialized by using the function twidfftrad2_fr16 for ifft_fr16 and twidfftrad2_fr32 for ifft_fr32. For optimal performance, the twiddle table should be allocated in a different memory section than the output array.
The argument twiddle_stride should be set to 1 if the twiddle table was originally created for an FFT of size fft_size. If the twiddle table was created for a larger FFT of size N*fft_size (where N is a power of 2), then twiddle_stride should be set to N. This argument therefore provides a way of using a single twiddle table to calculate FFTs of different sizes.
The argument scale_method controls how the function will apply scaling while computing a Fourier Transform. The available options are static scaling (dividing the input at any stage by 2), dynamic scaling (dividing the input at any stage by 2 if the largest absolute input value is greater or equal to 0.25), or no scaling. Note that the number of stages required to compute an FFT is dependent on the size of the FFT and is given by the formula log2(fft_size).
If static scaling is selected, the function will always scale intermediate results, thus preventing overflow. The loss of precision increases in line with fft_size and is more pronounced for input signals with a small magnitude (since the output is scaled by 1/fft_size). To select static scaling, set the argument scale_method to a value of 1. The block exponent returned will be log2(fft_size).
If dynamic scaling is selected, the function will inspect intermediate results and only apply scaling where required to prevent overflow. The loss of precision increases in line with the size of the FFT and is more pronounced for input signals with a large magnitude (since these factors increase the need for scaling). The requirement to inspect intermediate results will have an impact on performance. To select dynamic scaling, set the argument scale_method to a value of 2. The block exponent returned will be between 0 and log2(fft_size) depending upon the number of times that the function scales each set of intermediate results.
If no scaling is selected, the function will never scale intermediate results. There will be no loss of precision unless overflow occurs and in this case the function will generate saturated results. The likelihood of saturation increases in line with the FFT size and is more pronounced for input signals with a large magnitude. To select no scaling, set the argument scale_method to 3. The block exponent returned will be 0.
Note: Any values for the argument scale_method other than 2 or 3 will result in the function performing static scaling.
Error Conditions
The ifft functions abort if the FFT size is less than 8 or if the twiddle stride is less than 1.
Algorithm
The following equation is the basis of the algorithm.
N-1 -n*k
x(n)=Сумма X(k)*W
k=0 N
Domain
Input sequence length fft_size must be a power of 2 and at least 8.
Example
/* Compute IFFT( CFFT( X ) ) = X */
#include <filter.h>
#define N_FFT 64
complex_fract16 in[N_FFT];
complex_fract16 out_cfft[N_FFT];
complex_fract16 out_ifft[N_FFT];
complex_fract16 twiddle[N_FFT/2];
int blk_exp;
void ifft_fr16_example(void)
{
int i;
/* Generate DC signal */
for( i = 0; i < N_FFT; i++ )
{
in.re = 0x100;
in.im = 0x0;
}
/* Populate twiddle table */
twidfftrad2_fr16(twiddle, N_FFT);
/* Compute Fast Fourier Transform */
cfft_fr16(in, out_cfft, twiddle, 1, N_FFT, &blk_exp, 0);
/* Reverse static scaling applied by cfft_fr16() function
Apply the shift operation before the call to the
ifft_fr16() function only if all the values in out_cfft
= 0x100. Otherwise, perform the shift operation after the
ifft_fr16() function has been computed.
*/
for( i = 0; i < N_FFT; i++ )
{
out_cfft.re = out_cfft.re << 6; /* log2(N_FFT) = 6 */
out_cfft.im = out_cfft.im << 6;
}
/* Compute Inverse Fast Fourier Transform
The output signal from the ifft function will be the same
as the DC signal of magnitude 0x100 which was passed into
the cfft function.
*/
ifft_fr16(out_cfft, out_ifft, twiddle, 1, N_FFT, &blk_exp, 0);
}
N-point radix-2 inverse FFT
Synopsis
#include <filter.h>
void ifft_fr16(const complex_fract16 input[],
complex_fract16 output[],
const complex_fract16 twiddle_table[],
int twiddle_stride,
int fft_size,
int *block_exponent,
int scale_method);
void ifft_fr32(const complex_fract32 input[],
complex_fract32 output[],
const complex_fract32 twiddle_table[],
int twiddle_size,
int fft_size,
int *block_exponent,
int scale_method);
Description
The ifft functions transform the frequency domain complex input signal sequence to the time domain by using the radix-2 Fast Fourier Transform.
The size of the input array input and the output array output is fft_size, where fft_size represents the number of points in the FFT. By allocating these arrays in different memory banks, any potential data bank collisions are avoided, thus improving run-time performance. If the input data can be overwritten, the optimum memory usage can be achieved by also specifying the input array as the output array.
The twiddle table is passed in the argument twiddle_table, which must contain at least fft_size/2 twiddle factors. The table is composed of +cosine and -sine coefficients and may be initialized by using the function twidfftrad2_fr16 for ifft_fr16 and twidfftrad2_fr32 for ifft_fr32. For optimal performance, the twiddle table should be allocated in a different memory section than the output array.
The argument twiddle_stride should be set to 1 if the twiddle table was originally created for an FFT of size fft_size. If the twiddle table was created for a larger FFT of size N*fft_size (where N is a power of 2), then twiddle_stride should be set to N. This argument therefore provides a way of using a single twiddle table to calculate FFTs of different sizes.
The argument scale_method controls how the function will apply scaling while computing a Fourier Transform. The available options are static scaling (dividing the input at any stage by 2), dynamic scaling (dividing the input at any stage by 2 if the largest absolute input value is greater or equal to 0.25), or no scaling. Note that the number of stages required to compute an FFT is dependent on the size of the FFT and is given by the formula log2(fft_size).
If static scaling is selected, the function will always scale intermediate results, thus preventing overflow. The loss of precision increases in line with fft_size and is more pronounced for input signals with a small magnitude (since the output is scaled by 1/fft_size). To select static scaling, set the argument scale_method to a value of 1. The block exponent returned will be log2(fft_size).
If dynamic scaling is selected, the function will inspect intermediate results and only apply scaling where required to prevent overflow. The loss of precision increases in line with the size of the FFT and is more pronounced for input signals with a large magnitude (since these factors increase the need for scaling). The requirement to inspect intermediate results will have an impact on performance. To select dynamic scaling, set the argument scale_method to a value of 2. The block exponent returned will be between 0 and log2(fft_size) depending upon the number of times that the function scales each set of intermediate results.
If no scaling is selected, the function will never scale intermediate results. There will be no loss of precision unless overflow occurs and in this case the function will generate saturated results. The likelihood of saturation increases in line with the FFT size and is more pronounced for input signals with a large magnitude. To select no scaling, set the argument scale_method to 3. The block exponent returned will be 0.
Note: Any values for the argument scale_method other than 2 or 3 will result in the function performing static scaling.
Error Conditions
The ifft functions abort if the FFT size is less than 8 or if the twiddle stride is less than 1.
Algorithm
The following equation is the basis of the algorithm.
N-1 -n*k
x(n)=Сумма X(k)*W
k=0 N
Domain
Input sequence length fft_size must be a power of 2 and at least 8.
Example
/* Compute IFFT( CFFT( X ) ) = X */
#include <filter.h>
#define N_FFT 64
complex_fract16 in[N_FFT];
complex_fract16 out_cfft[N_FFT];
complex_fract16 out_ifft[N_FFT];
complex_fract16 twiddle[N_FFT/2];
int blk_exp;
void ifft_fr16_example(void)
{
int i;
/* Generate DC signal */
for( i = 0; i < N_FFT; i++ )
{
in.re = 0x100;
in.im = 0x0;
}
/* Populate twiddle table */
twidfftrad2_fr16(twiddle, N_FFT);
/* Compute Fast Fourier Transform */
cfft_fr16(in, out_cfft, twiddle, 1, N_FFT, &blk_exp, 0);
/* Reverse static scaling applied by cfft_fr16() function
Apply the shift operation before the call to the
ifft_fr16() function only if all the values in out_cfft
= 0x100. Otherwise, perform the shift operation after the
ifft_fr16() function has been computed.
*/
for( i = 0; i < N_FFT; i++ )
{
out_cfft.re = out_cfft.re << 6; /* log2(N_FFT) = 6 */
out_cfft.im = out_cfft.im << 6;
}
/* Compute Inverse Fast Fourier Transform
The output signal from the ifft function will be the same
as the DC signal of magnitude 0x100 which was passed into
the cfft function.
*/
ifft_fr16(out_cfft, out_ifft, twiddle, 1, N_FFT, &blk_exp, 0);
}
-
Бахурин Сергей
- Администратор
- Сообщения: 1116
- Зарегистрирован: 05 окт 2010, 19:55
- Контактная информация:
Re: БПФ
давайте разберемся.
Прямое преобразование:
)
Когда вы делали fixedpoint прямое fft вы сдвигами вправо выравнивали положение точки правильно я понимаю? Т.е. на выходе прямого преобразования вы получили правильные значения спектра.
Обратное преобразование:
)
Для выравнивания положения точки обратного преобразования вы также использовали сдвиги вправо. Если вы все правильно сделали то результат обратного преобразования необходимо поделить на N (которая в формуле для обратного преобразования перед суммой) т.е. сделать 9 сдвигов вправо один раз поскольку положение точки верно в обоих случаях. Если я прав то без деления на N на выходе ifft сигнал будет в 512 раз выше исходного.
Прямое преобразование:
Когда вы делали fixedpoint прямое fft вы сдвигами вправо выравнивали положение точки правильно я понимаю? Т.е. на выходе прямого преобразования вы получили правильные значения спектра.
Обратное преобразование:
Для выравнивания положения точки обратного преобразования вы также использовали сдвиги вправо. Если вы все правильно сделали то результат обратного преобразования необходимо поделить на N (которая в формуле для обратного преобразования перед суммой) т.е. сделать 9 сдвигов вправо один раз поскольку положение точки верно в обоих случаях. Если я прав то без деления на N на выходе ifft сигнал будет в 512 раз выше исходного.
Re: БПФ
По теории все верно. Но у меня после обратного преобразования сигнал получается в 512 раз меньше исходного. Это с одной стороны. Теперь относительно реализации для BlackFin - выше я привел пример их кода(из Help-а для функции прямого преобразования) и там идет сравнение с числом 0x100 для 16 разрядных чисел с фиксированной точкой. Там же приводится рекомендация по масштабированию сигнала после прямого преобразования:
/* Reverse static scaling applied by cfft_fr16() function
Apply the shift operation before the call to the
ifft_fr16() function only if all the values in out_cfft
= 0x100. Otherwise, perform the shift operation after the
ifft_fr16() function has been computed.
*/
и далее
/* Compute Inverse Fast Fourier Transform
The output signal from the ifft function will be the same
as the DC signal of magnitude 0x100 which was passed into
the cfft function.
*/
непонятно почему так...
/* Reverse static scaling applied by cfft_fr16() function
Apply the shift operation before the call to the
ifft_fr16() function only if all the values in out_cfft
= 0x100. Otherwise, perform the shift operation after the
ifft_fr16() function has been computed.
*/
и далее
/* Compute Inverse Fast Fourier Transform
The output signal from the ifft function will be the same
as the DC signal of magnitude 0x100 which was passed into
the cfft function.
*/
непонятно почему так...
-
Бахурин Сергей
- Администратор
- Сообщения: 1116
- Зарегистрирован: 05 окт 2010, 19:55
- Контактная информация:
Re: БПФ
Я думаю надо повнимательнее поглядеть что лежит в таблицах синусов и косинусов а именно в twiddle после вызова
Поскольку синусы и косинусы меньше 1 а математика в интах то надо понять где учитывается положение точки при преобразовании.
Еще вопрос такой что выдает функция прямого преобразования в массиве out_cfft перед этим циклом?
с заданным входным сигналом 0x100 и N_FFT = 64? От этого можно отплясать чтобы понять как правильно нормировать результат прямого и обратного преобразования. Должно быть out_cfft[0].re = 0x4000 все остальное нули.
Код: Выделить всё
/* Populate twiddle table */
twidfftrad2_fr16(twiddle, N_FFT);
Еще вопрос такой что выдает функция прямого преобразования в массиве out_cfft перед этим циклом?
Код: Выделить всё
for( i = 0; i < N_FFT; i++ )
{
out_cfft[i].re = out_cfft[i].re << 6; /* log2(N_FFT) = 6 */
out_cfft[i].im = out_cfft[i].im << 6;
}
Re: БПФ
В массиве out_cfft функция cfft_fr16 выдант frequency domain complex output signal sequence (последовательность в частотной области в виде комплексного массива, т.е спектр сигнала).
Сигнал 0x100 это прямоугольный импульс длительностью 64 точки. А вот почему они пишут
Apply the shift operation before the call to the ifft_fr16() function only if all the values in out_cfft = 0x100 (выполните рперацию сдвига перед вызовом обратного преобразования если все значения в out_cfft = 0x100) - это половина разрядной сетки 16 разрядного числа ??? Получается они проверяют перед масштабированием не произойдет ли переполнение...
Сигнал 0x100 это прямоугольный импульс длительностью 64 точки. А вот почему они пишут
Apply the shift operation before the call to the ifft_fr16() function only if all the values in out_cfft = 0x100 (выполните рперацию сдвига перед вызовом обратного преобразования если все значения в out_cfft = 0x100) - это половина разрядной сетки 16 разрядного числа ??? Получается они проверяют перед масштабированием не произойдет ли переполнение...
-
Бахурин Сергей
- Администратор
- Сообщения: 1116
- Зарегистрирован: 05 окт 2010, 19:55
- Контактная информация:
Re: БПФ
Что лежит в массиве out_cfft перед циклом?
Re: БПФ
спектр сигнала...гармоники и фазы
Спектр аналитического сигнала:
еще вопрос по преобразованию Гилберта(я его задал в теме по преобразованию Гилберта, но вы не ответили-кто то дал ответ который я не просил)
Исходный сигнал подвергнуть преобразованию Фурье, обнулить спектр в отрицательной области частот, удвоить спектр в положительной области частот, после взять обратное преобразование Фурье и получится аналитический сигнал, из которого можно выделить исходный сигнал и его ортогональное дополнение.
А положительную частоту с индексом 0 ( S(0), w=0 ) не надо удваивать ?
Спектр аналитического сигнала:
еще вопрос по преобразованию Гилберта(я его задал в теме по преобразованию Гилберта, но вы не ответили-кто то дал ответ который я не просил)
Исходный сигнал подвергнуть преобразованию Фурье, обнулить спектр в отрицательной области частот, удвоить спектр в положительной области частот, после взять обратное преобразование Фурье и получится аналитический сигнал, из которого можно выделить исходный сигнал и его ортогональное дополнение.
А положительную частоту с индексом 0 ( S(0), w=0 ) не надо удваивать ?
-
Бахурин Сергей
- Администратор
- Сообщения: 1116
- Зарегистрирован: 05 окт 2010, 19:55
- Контактная информация:
Re: БПФ
1. Я понимаю что спектр что именно лежит в массиве. какие именно числа?
2. Преобразование Гильберта для нулевой частоты не надо удваивать. А лучше выбросить постоянку из сигнала когда делаете преобразование Гильберта.
2. Преобразование Гильберта для нулевой частоты не надо удваивать. А лучше выбросить постоянку из сигнала когда делаете преобразование Гильберта.
Re: БПФ
Да, все верно, в массиве out_cfft[0].re = 0x4000 все остальное нули, в out_cfft.im все нули.