DFT in C ...

Тема в разделе "WASM.AUDIO", создана пользователем Andrei, 3 май 2018.

  1. Andrei

    Andrei Member

    Публикаций:
    0
    Регистрация:
    13 апр 2018
    Сообщения:
    322
    Почему уровень сигнала на спектре зависит от частоты? При одинаковом сигнале на входе ???
     
  2. Thetrik

    Thetrik UA6527P

    Публикаций:
    0
    Регистрация:
    25 июл 2011
    Сообщения:
    868
    Если частота кратна окну то не зависит.
     
  3. Andrei

    Andrei Member

    Публикаций:
    0
    Регистрация:
    13 апр 2018
    Сообщения:
    322
    "ОКНО" - Это длинна буфера который участвует в вычислениях ?
    Я правильно понимаю ?
     
  4. Andrei

    Andrei Member

    Публикаций:
    0
    Регистрация:
    13 апр 2018
    Сообщения:
    322
    Вообще в принципе в коде обработки DFT нет частотозадающих элементов.
    Я так понял что частота определяется нарезкой входных семплов, от АЦП ...
     
  5. Andrei

    Andrei Member

    Публикаций:
    0
    Регистрация:
    13 апр 2018
    Сообщения:
    322
    Ты как профи вот скажи вот этот код можно использовать ????

    Код (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.     }
     
    Последнее редактирование модератором: 3 май 2018
  6. SadKo

    SadKo Владимир Садовников

    Публикаций:
    8
    Регистрация:
    4 июн 2007
    Сообщения:
    1.610
    Адрес:
    г. Санкт-Петербург
    rand(), seriously?
     
  7. Thetrik

    Thetrik UA6527P

    Публикаций:
    0
    Регистрация:
    25 июл 2011
    Сообщения:
    868
    Преобразование Фурье работает только с бесконечными сигналами. К примеру если берется прямоугольное окно (нет окна), то это эквивалентно свертке в частотной области с частотной характеристикой окна (sinc). Если к примеру у нас чистая синусоида во временной области, то это соответствует дельта-функции в частотной, а свертка с дельта-функцией эквивалентно сдвигу. Если период синусоиды точно кратен размеру окна, то "хвосты" ликвидируются, если синусоида расположена где-то между частотными выборками соответственно эти "хвосты накладываются". Накидал небольшую анимацию с объяснением принципа:
    [​IMG]
    Поэтому если синусоида не "попадает в сетку", то она размазывается из-за этого амплитуда зависит от частоты. При увеличении размера окна, оконная функция стремится к дельта-функции. Для уменьшения амплитуды хвостов используются различные весовые окна.
     
    Andrei нравится это.
  8. Thetrik

    Thetrik UA6527P

    Публикаций:
    0
    Регистрация:
    25 июл 2011
    Сообщения:
    868
    Частота в дискретном преобразовании Фурье задается в долях от частоты дискретизации, от 0 до 0.5D, после половины частоты дискретизации спектр реального дискретного сигнала зеркально симметричен и бесконечно периодичен в обе стороны.

    Для чего?
     
    Andrei нравится это.
  9. SadKo

    SadKo Владимир Садовников

    Публикаций:
    8
    Регистрация:
    4 июн 2007
    Сообщения:
    1.610
    Адрес:
    г. Санкт-Петербург
    Для генерации белого шума, видимо.
     
  10. Andrei

    Andrei Member

    Публикаций:
    0
    Регистрация:
    13 апр 2018
    Сообщения:
    322
    Хочу сделать примитивный анализатор спектра, ...
     
  11. Andrei

    Andrei Member

    Публикаций:
    0
    Регистрация:
    13 апр 2018
    Сообщения:
    322
    Генератор белого шума+sin() это чтобы смоделировать входной сигнал на фоне шума )))
     
  12. Andrei

    Andrei Member

    Публикаций:
    0
    Регистрация:
    13 апр 2018
    Сообщения:
    322
    Насчет уровня понял , насчет окна не совсем, вот в этом коде где ОКНО ? И где задаётся его размер ??

    Код (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.    }
     
    Последнее редактирование модератором: 4 май 2018
  13. Andrei

    Andrei Member

    Публикаций:
    0
    Регистрация:
    13 апр 2018
    Сообщения:
    322
    Еще не понятно вот что на входе уровень в Вольтах, уровень на выходе тоже в Вольтах или Ваттах ??? или в чем ?
     
  14. SadKo

    SadKo Владимир Садовников

    Публикаций:
    8
    Регистрация:
    4 июн 2007
    Сообщения:
    1.610
    Адрес:
    г. Санкт-Петербург
    Если в вашем коде не применяется оконная функция, то считайте, что используется прямоугольная оконная функция.
    Поподробнее можно прочитать здесь:
    https://en.wikipedia.org/wiki/Window_function
     
  15. Andrei

    Andrei Member

    Публикаций:
    0
    Регистрация:
    13 апр 2018
    Сообщения:
    322
    Еще интересует как правильно сделать эквалайзер на три канала НЧ,ВЧ,СЧ ...

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

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

    Andrei Member

    Публикаций:
    0
    Регистрация:
    13 апр 2018
    Сообщения:
    322
    Так "Окно" - ЭТО интервал между соседними отсчетами идущими от АЦП ???
     
  17. Thetrik

    Thetrik UA6527P

    Публикаций:
    0
    Регистрация:
    25 июл 2011
    Сообщения:
    868
    Нет. Я же написал ПФ берется от бесконечного сигнала и также преобразует в бесконечный спектр. Когда ты берешь только часть данных, то считай что ты берешь бесконечный периодический сигнал и умножаешь его на прямоугольник, остальные данные зануляешь.
     
  18. Thetrik

    Thetrik UA6527P

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

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

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

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

    Thetrik UA6527P

    Публикаций:
    0
    Регистрация:
    25 июл 2011
    Сообщения:
    868
    https://wasm.in/threads/gotovye-proekty-na-vb6.31728/#post-383200

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

    Andrei Member

    Публикаций:
    0
    Регистрация:
    13 апр 2018
    Сообщения:
    322
    Книга по DSP действительно замечательная, ... мне до пенсии хватит Ёе читать ))))
     
    Последнее редактирование: 5 май 2018