DFT in C ...

Почему уровень сигнала на спектре зависит от частоты? При одинаковом сигнале на входе ???

 

Если частота кратна окну то не зависит.

 

"ОКНО" - Это длинна буфера который участвует в вычислениях ?
Я правильно понимаю ?

 

Вообще в принципе в коде обработки DFT нет частотозадающих элементов.
Я так понял что частота определяется нарезкой входных семплов, от АЦП ...

 

Ты как профи вот скажи вот этот код можно использовать ????

Код (C++):

1. #include <stdio.h>
2. #include <stdlib.h>
3. #include <math.h>
4. // N is assumed to be greater than 4 and a power of 2
5. #define N 64
6. #define PI2 6.2832
7. // Twiddle factors (64th roots of unity)
8. const float W[] = {
9.      1.00000, 0.99518, 0.98079, 0.95694, 0.92388, 0.88192, 0.83147, 0.77301,
10.      0.70711, 0.63439, 0.55557, 0.47139, 0.38268, 0.29028, 0.19509, 0.09801,
11.     -0.00000,-0.09802,-0.19509,-0.29029,-0.38269,-0.47140,-0.55557,-0.63440,
12.     -0.70711,-0.77301,-0.83147,-0.88192,-0.92388,-0.95694,-0.98079,-0.99519,
13.     -1.00000,-0.99518,-0.98078,-0.95694,-0.92388,-0.88192,-0.83146,-0.77300,
14.     -0.70710,-0.63439,-0.55556,-0.47139,-0.38267,-0.29027,-0.19508,-0.09801,
15.      0.00001, 0.09803, 0.19510, 0.29030, 0.38269, 0.47141, 0.55558, 0.63440,
16.      0.70712, 0.77302, 0.83148, 0.88193, 0.92388, 0.95694, 0.98079, 0.99519
17. };
18. int main()
19. {
20.     // time and frequency domain data arrays
21.     int n, k;                     // time and frequency domain indices
22.     float x[N];                   // discrete-time signal, x
23.     float Xre[N/2+1], Xim[N/2+1]; // DFT of x (real and imaginary parts)
24.     float P[N/2+1];               // power spectrum of x
25.     // Generate random discrete-time signal x in range (-1,+1)
26.     srand(time(0));
27.     for (n=0 ; n<N ; ++n) x[n] = (((2.0 * rand()) / RAND_MAX) - 1.0+sin(PI2*n*3.7/N))/2;
28.    // for (n=0 ; n<N ; ++n) x[n] = ((2.0 * rand()) / RAND_MAX) - 1.0;
29.     // Calculate DFT and power spectrum up to Nyquist frequency
30.     int to_sin = 3*N/4; // index offset for sin
31.     int a, b;
32.     for (k=0 ; k<=N/2 ; ++k)
33.     {
34.         Xre[k] = 0; Xim[k] = 0;
35.         a = 0; b = to_sin;
36.         for (n=0 ; n<N ; ++n)
37.         {
38.             Xre[k] += x[n] * W[a%N];
39.             Xim[k] -= x[n] * W[b%N];
40.             a += k; b += k;
41.         }
42.         P[k] = Xre[k]*Xre[k] + Xim[k]*Xim[k];
43.     }

 

rand(), seriously?

 

Преобразование Фурье работает только с бесконечными сигналами. К примеру если берется прямоугольное окно (нет окна), то это эквивалентно свертке в частотной области с частотной характеристикой окна (sinc). Если к примеру у нас чистая синусоида во временной области, то это соответствует дельта-функции в частотной, а свертка с дельта-функцией эквивалентно сдвигу. Если период синусоиды точно кратен размеру окна, то "хвосты" ликвидируются, если синусоида расположена где-то между частотными выборками соответственно эти "хвосты накладываются". Накидал небольшую анимацию с объяснением принципа:

Поэтому если синусоида не "попадает в сетку", то она размазывается из-за этого амплитуда зависит от частоты. При увеличении размера окна, оконная функция стремится к дельта-функции. Для уменьшения амплитуды хвостов используются различные весовые окна.

 

Частота в дискретном преобразовании Фурье задается в долях от частоты дискретизации, от 0 до 0.5D, после половины частоты дискретизации спектр реального дискретного сигнала зеркально симметричен и бесконечно периодичен в обе стороны.

Для чего?

 

Для генерации белого шума, видимо.

 

Хочу сделать примитивный анализатор спектра, ...

 

Генератор белого шума+sin() это чтобы смоделировать входной сигнал на фоне шума )))

 

Насчет уровня понял , насчет окна не совсем, вот в этом коде где ОКНО ? И где задаётся его размер ??

Код (C++):

1. #include <stdio.h>
2. #include <stdlib.h>
3. #include <math.h>
4. // N is assumed to be greater than 4 and a power of 2
5. #define N 64
6. #define PI2 6.2832
7. // Twiddle factors (64th roots of unity)
8. const float W[] = {
9.      1.00000, 0.99518, 0.98079, 0.95694, 0.92388, 0.88192, 0.83147, 0.77301,
10.      0.70711, 0.63439, 0.55557, 0.47139, 0.38268, 0.29028, 0.19509, 0.09801,
11.    -0.00000,-0.09802,-0.19509,-0.29029,-0.38269,-0.47140,-0.55557,-0.63440,
12.    -0.70711,-0.77301,-0.83147,-0.88192,-0.92388,-0.95694,-0.98079,-0.99519,
13.    -1.00000,-0.99518,-0.98078,-0.95694,-0.92388,-0.88192,-0.83146,-0.77300,
14.    -0.70710,-0.63439,-0.55556,-0.47139,-0.38267,-0.29027,-0.19508,-0.09801,
15.      0.00001, 0.09803, 0.19510, 0.29030, 0.38269, 0.47141, 0.55558, 0.63440,
16.      0.70712, 0.77302, 0.83148, 0.88193, 0.92388, 0.95694, 0.98079, 0.99519
17. };
18. int main()
19. {
20.    // time and frequency domain data arrays
21.    int n, k;                   // time and frequency domain indices
22.    float x[N];                 // discrete-time signal, x
23.    float Xre[N/2+1], Xim[N/2+1]; // DFT of x (real and imaginary parts)
24.    float P[N/2+1];             // power spectrum of x
25.    // Generate random discrete-time signal x in range (-1,+1)
26.    srand(time(0));
27.    for (n=0 ; n<N ; ++n) x[n] = (((2.0 * rand()) / RAND_MAX) - 1.0+sin(PI2*n*3.7/N))/2;
28.    // for (n=0 ; n<N ; ++n) x[n] = ((2.0 * rand()) / RAND_MAX) - 1.0;
29.    // Calculate DFT and power spectrum up to Nyquist frequency
30.    int to_sin = 3*N/4; // index offset for sin
31.    int a, b;
32.    for (k=0 ; k<=N/2 ; ++k)
33.    {
34.         Xre[k] = 0; Xim[k] = 0;
35.         a = 0; b = to_sin;
36.        for (n=0 ; n<N ; ++n)
37.        {
38.             Xre[k] += x[n] * W[a%N];
39.             Xim[k] -= x[n] * W[b%N];
40.             a += k; b += k;
41.        }
42.         P[k] = Xre[k]*Xre[k] + Xim[k]*Xim[k];
43.    }

 

Еще не понятно вот что на входе уровень в Вольтах, уровень на выходе тоже в Вольтах или Ваттах ??? или в чем ?

 

Если в вашем коде не применяется оконная функция, то считайте, что используется прямоугольная оконная функция.
Поподробнее можно прочитать здесь:
https://en.wikipedia.org/wiki/Window_function

 

Еще интересует как правильно сделать эквалайзер на три канала НЧ,ВЧ,СЧ ...

1. Берем три фильтра
2. В каждом фильтре своё усиление
3. Потом складываем

Фильтры какого порядка нужны ?
Как сложить (кан1+кан2+кан3)/3

 

Так "Окно" - ЭТО интервал между соседними отсчетами идущими от АЦП ???

 

Нет. Я же написал ПФ берется от бесконечного сигнала и также преобразует в бесконечный спектр. Когда ты берешь только часть данных, то считай что ты берешь бесконечный периодический сигнал и умножаешь его на прямоугольник, остальные данные зануляешь.

 

Можешь реализовать готовый фильтр уже с необходимой АЧХ. Самый простой вариант, нарисовать АЧХ (ФЧХ), выполнить ОПФ, полученную импульсную характеристику умножить на сглаживающее окно (если после этого выполнить ПФ, то узнаешь реальную АЧХ и ФЧХ) и использовать ее в качестве ядра КИХ фильтра.

Если качество (АЧХ) не сильно важно то можно использовать MA фильтр, у него очень эффективная рекурсивная реализация, но плохая АЧХ (sinc), если применить фильтр несколько раз, то можно улучшить АЧХ (в этом случае АЧХ умножаются, а ИХ свертываются, получается прямоугольник, треугольник, и что-то напоминающее гауссиан). Благодаря тому что MA фильтр имеет линейную ФЧХ то другие типы фильтров реализуются простым вычитанием сигналов на этапах фильтрации между собой.

Также можно рассчитать БИХ фильтр на основе нулей и полюсов, а частотную характеристику потом получить из z-преобразования на единичной окружности, но т.к. этот тип фильтров нестабилен, то проще пропустить дельта импульс через этот фильтр и посмотреть его АЧХ.

Вот готовое описание эквалайзера на all-pass фильтре. Вот тут куча готовых реализаций. Замечательная книга по DSP, все доступным языком описано.

 

https://wasm.in/threads/gotovye-proekty-na-vb6.31728/#post-383200

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

 

Книга по DSP действительно замечательная, ... мне до пенсии хватит Ёе читать ))))