Форум программистов «Весельчак У»
  *
Добро пожаловать, Гость. Пожалуйста, войдите или зарегистрируйтесь.
Вам не пришло письмо с кодом активации?

  • Рекомендуем проверить настройки временной зоны в вашем профиле (страница "Внешний вид форума", пункт "Часовой пояс:").
  • У нас больше нет рассылок. Если вам приходят письма от наших бывших рассылок mail.ru и subscribe.ru, то знайте, что это не мы рассылаем.
   Начало  
Наши сайты
Помощь Поиск Календарь Почта Войти Регистрация  
 
Страниц: [1]   Вниз
  Печать  
Автор Тема: Быстрое преобразование Фурье  (Прочитано 15849 раз)
0 Пользователей и 1 Гость смотрят эту тему.
sinsin
Постоялец

ru
Offline Offline
Пол: Мужской

« : 24-02-2011 10:31 » 

Доброго времени суток.
Пытаюсь реализовать на С++ создание цифрового фильтра методом свертки, при этом использую БПФ, реализованную на http://algolist.manual.ru/maths/fft.php, но получаю странный результат:
после прогона заданной частотной характеристики фильтра через обратное БПФ(пока использую FFT_T) в выходном массиве появляеются ненулевые комплексные составляющие.
Код:
           int N = FDiscr*2;
Complex *FHar = new Complex[N];
double *tmp_arr = new double[N];
Complex tmp;
// частотная характеристика
for (int i = 0;i<FDiscr;i++)
{
FHar[i].i = 0;
if ((i<=f2)&&(i>=f1))
FHar[i].r = 1;
else
FHar[i].r = 0;
}
FFTReOrder(FHar,FDiscr);
// обратное бпф
FFT_T(FHar,FDiscr,-1);

Как это можно побороть?
Чтобы насчет правильности действий не возникалли вопросы, синтез фильтра в mathcad прикреплен в Filter_exper.xmcd.

* Filter_exper.xmcd (122.65 Кб - загружено 1097 раз.)
« Последнее редактирование: 24-02-2011 11:13 от sinsin » Записан
RXL
Технический
Администратор

ru
Offline Offline
Пол: Мужской

WWW
« Ответ #1 : 24-02-2011 12:47 » 

Не берусь разбираться в алгоритме, но хочу напомнить, что плавающие вычисления всегда будут содержать погрешности - промежуточный и/или конечный результат надо округлять до требуемой точности.
Записан

... мы преодолеваем эту трудность без синтеза распределенных прототипов. (с) Жуков М.С.
Dimka
Деятель
Команда клуба

ru
Offline Offline
Пол: Мужской

« Ответ #2 : 24-02-2011 13:07 » 

sinsin, советую решать задачу по частям. В начале отладь БПФ и обратное БПФ, потом уже занимайся частотным фильтром. Для этого возьми массив комплексных значений (можно с нулевыми мнимыми частями), выполни для него БПФ и получи спектр, потом для спектра выполни и обратное БПФ и добивайся получения исходного массива (с разумной погрешностью).
Записан

Программировать - значит понимать (К. Нюгард)
Невывернутое лучше, чем вправленное (М. Аврелий)
Многие готовы скорее умереть, чем подумать (Б. Рассел)
sinsin
Постоялец

ru
Offline Offline
Пол: Мужской

« Ответ #3 : 24-02-2011 13:58 » 

RXL, там слишком большие погрешности получаются. Dimka, библиотека вроде рабочая, мне её порекомендовали как уже использующуюся. Завтра еще покопаю,может что найду.
Записан
Ochkarik
Команда клуба

ru
Offline Offline
Пол: Мужской

« Ответ #4 : 24-02-2011 17:15 » 

sinsin, а вы как хотели? у вас фильтр - не симметричный в частотной области относительно 0 частоты.
 само собой - у него комплексная ИХ будет.
FDiscr - это размер БПФ или это половина? ничего не понял.
маткада под рукой нет, посмотреть не могу...
Записан

RTFM уже хоть раз наконец!  RTFM :[ ну или хотя бы STFW...
sinsin
Постоялец

ru
Offline Offline
Пол: Мужской

« Ответ #5 : 24-02-2011 17:22 » 

Ochkarik, размер.
Записан
Ochkarik
Команда клуба

ru
Offline Offline
Пол: Мужской

« Ответ #6 : 24-02-2011 17:26 » 

чтобы получить чисто действительную ИХ фильтра - частотная характеристика фильтра должна быть симметричной относительно нулевой частоты (или Fд/2, что то же самое).
так что просто сделайте зеркальное отображение

Добавлено через 2 минуты и 21 секунду:
PS зачем тогда у вас массивы удвоенные FDiscr*2?
« Последнее редактирование: 24-02-2011 17:28 от Ochkarik » Записан

RTFM уже хоть раз наконец!  RTFM :[ ну или хотя бы STFW...
sinsin
Постоялец

ru
Offline Offline
Пол: Мужской

« Ответ #7 : 24-02-2011 17:45 » 

Ochkarik, спасибо за совет. Для резервирования: дальше по методичке (http://graphicon.ru/oldgr/courses/cg02b/library/dspcourse.pdf, стр 27-28) отражаю массив относительно ноля (можно сделать и "на месте", но пока не нужно, хочу просто фильтр отладить), отбрасываю отсчеты, применяю функцию окна (Хэмминг)

Добавлено через 12 часов, 58 минут и 20 секунд:
Ochkarik, такой вопрос: частотная характеристика в Mathcad'е задана точно также, как и в коде, но получаю действительную импульсную характеристику, вчем может быть подвох? (пересохранил маткадовский файл в rtf)

Добавлено через 23 часа, 49 минут и 30 секунд:
Нашел свою ошибку, частотную характеристику фильтра надо задачать так:
Код:
              FHar[0].m_im = FHar[FDiscr].m_im = 0;
FHar[0].m_re = (f1==0)?1:0;
FHar[FDiscr].m_re = (f2==FDiscr)?1:0;
// частотная характеристика
for (int i = 1;i<FDiscr;i++)
{
FHar[i].m_im = 0;
if ((i<=f2)&&(i>=f1))
FHar[i].m_re = 1;
else
FHar[i].m_re = 0;
FHar[N-i] = FHar[i];
}
Все фильтрует отлично, однако возникают проблемы с нормированием сиргнала: при классическом делении каждой точки  обратного БПФ амплитуда синусоиды с единичной становится 0,077. Есть ли какие-нибудь еще способы нормирования?

* Filter_exper.rtf (377.27 Кб - загружено 1113 раз.)
« Последнее редактирование: 26-02-2011 06:33 от sinsin » Записан
RXL
Технический
Администратор

ru
Offline Offline
Пол: Мужской

WWW
« Ответ #8 : 26-02-2011 09:18 » 

sinsin, покажи пример - что-то непонятно, как 1 в 0.077 преобразуется.
Записан

... мы преодолеваем эту трудность без синтеза распределенных прототипов. (с) Жуков М.С.
sinsin
Постоялец

ru
Offline Offline
Пол: Мужской

« Ответ #9 : 26-02-2011 09:32 » 

Формирование массива происходит след. образом:
Код:
             for (int i = 0;i<65536;i++)
{
In[i] = 2+sin(120*2*M_PI*i/N)+5*sin(2*M_PI*i/N);
}
Filter(In, Out, 0, 100, 65536/2);
В функции Filter параметры: In - вх. сигнал, Out - выход, след. два параметра - частоты среза фильтра (в данном случае имеем ФНЧ, "отрезающий" все, выше 100 герц).
Как я понимаю, я должен на выходе получить сигнал в виде суммы 2+5*sin(2*M_PI*i/N).
Выход фильтра показан на рисунке.

Т.е фильтрация есть, но она не нормирована.
Сейчас нормирую в двух местах:
1) при обратном БПФ  делю все отсчеты на N;
2) делю каждый отсчет импульсной характеристики фильтра после наложения окна на сумму модулей её составляющих (там получаются действительные числа, мнимая часть порядка 10^-4 - 10^-5 - погрешность).
Почему думаю на нормализацию: множитель получается не константным, в частности выше он где-то около 5, а для  массива, сформированного
Код:
             for (int i = 0;i<65536;i++)
{
In[i] = 2+sin(120*2*M_PI*i/N)+5*sin(20*2*M_PI*i/N);
}
Filter(In, Out, 10, 100, 65536/2);
он около 10

* Filter.PNG (38.43 Кб - загружено 1045 раз.)
* Filter.PNG (25.8 Кб - загружено 2536 раз.)
* Filter2.PNG (47.78 Кб - загружено 2491 раз.)
« Последнее редактирование: 26-02-2011 10:17 от sinsin » Записан
Ochkarik
Команда клуба

ru
Offline Offline
Пол: Мужской

« Ответ #10 : 26-02-2011 10:48 » 

2 - не есть правильно.  сейчас с ходу не соображу, но по моему коэф. передачи для наложенного окна надо считать отдельно от ИХ фильтра.

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

Добавлено через 9 минут и 44 секунды:
да, и коэффициент передачи считается не через сумму модулей. просто сумма(это будет коэф передачи на нулевой частоте). для расчета передачи оконной функции думаю подойдет. для фильтра в целом - нет.
 причем мнимая часть у вас должна получится нулевая, откуда там погрешность, если ИХ - чисто действительная?
или я что то не понял?
« Последнее редактирование: 26-02-2011 10:58 от Ochkarik » Записан

RTFM уже хоть раз наконец!  RTFM :[ ну или хотя бы STFW...
sinsin
Постоялец

ru
Offline Offline
Пол: Мужской

« Ответ #11 : 26-02-2011 12:35 » new

Ochkarik,  там действительно должны быть только действительные числа, погрешность  получается из-за плавающих вычислений, см #1

Добавлено через 1 час, 37 минут и 46 секунд:
Ochkarik, Вы правы, это из-за окна, если убрать окно, то получается все в порядке.
Нашел ресурс, где описан такой-же синтез фильтра, может кому-нибудь пригодится: http://dmilvdv.narod.ru/SpeechSynthesis/index.html?sound_info.html, тут, кстати , нормирование как было в пункте 2
« Последнее редактирование: 26-02-2011 14:12 от sinsin » Записан
Страниц: [1]   Вверх
  Печать  
 

Powered by SMF 1.1.21 | SMF © 2015, Simple Machines