Быстрое преобразование Фурье

[sinsin]

Доброго времени суток.
Пытаюсь реализовать на С++ создание цифрового фильтра методом свертки, при этом использую БПФ, реализованную на 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.

[RXL]

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

[Dimka]

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

[sinsin]

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

[Ochkarik]

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

[sinsin]

Ochkarik, размер.

[Ochkarik]

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

Добавлено через 2 минуты и 21 секунду:
PS зачем тогда у вас массивы удвоенные FDiscr*2?

[sinsin]

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. Есть ли какие-нибудь еще способы нормирования?

[RXL]

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

[sinsin]

Формирование массива происходит след. образом:

Код:

              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

[Ochkarik]

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

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

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

[sinsin]

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

Добавлено через 1 час, 37 минут и 46 секунд:
Ochkarik, Вы правы, это из-за окна, если убрать окно, то получается все в порядке.
Нашел ресурс, где описан такой-же синтез фильтра, может кому-нибудь пригодится: http://dmilvdv.narod.ru/SpeechSynthesis/index.html?sound_info.html, тут, кстати , нормирование как было в пункте 2