Устойчивость БПФ/ОБПФ на "длинных" сигналах

[Ivan Karamazov]

Про последний тезис -- более-менне в курсе. Речь только о возможном "развале" длинного сигнала из-за ошибок округления при ДПФ.
А "насчет зависимости от характера сигнала" -- все оказалось очень просто (как и предполагалось) -- в тестовом сигнале присутствовало одно локальное ограничение ("clip") на 4-5 остсчетов. Это дало "выбоину" на одном отсчете (причем только у мнимой части) Гильберта. Т.е. не фатальную выбоину -- просто на несколько младших разрядов у разных реализаций БПФ. (Ну да, избегайте ограничений -- они имеют "богатый" спектр со всеми вытекющими). И, как уже упоминалось -- если ваш сигнал имеет четное число отсчетов (при необходимости его легко разумным образом "доопределить") -- необходимо обнулить в т.ч. и n/2-ю компоненту спектра.

[Santik]

А сколько пар преобразований (FFT->IFFT) -> (FFT->IFFT) ->...можно сделать без "развала" исходного сигнала длиной N?
Можно оценить Nmax=N*2^M

[Ivan Karamazov]

Вообще, такой эксперимент дорог и нечист, хотя, конечно, попробовать можно. На самом деле я сделал так -- ввел режим верификации, выводящий
(1) статистику разности действительной части комплексного сигнала и исходного сигнала;
(2) статистику подавления отрицательных частот в спектре комплексного сигнала (т.е. для этого еще одно ДПФ делается).
Результаты -- разница в действительных частях в пределах младшего (из исходных 16-ти) бит, подавление -- не хуже -170 дБ (обычно порядка 200 (!) дБ. Но это -- дорого, поскольку у меня всего 4 ГБ RAM и Win7, у которой после однократного перерасхода памяти любит возникать ситуация: мой процесс в данный момент потребляет всего 2,7 ГБ, а размер дискового кэша -- 3,5 ГБ! Естественно, все стоит, а кэш общается с диском. (Такую систему угробили )

[Santik]

Вообще, такой эксперимент дорог и нечист

Извините, опять "по диагонали" все прочитал. Но чтобы Вашу методику осознать мне потребуется очень много времени.
"дорог"
Моя мысль в том , что необязательно брать N= 2^27. Взять 2^16 и много раз прокрутить. Памяти много не надо на это. И следить за мнимой частью выходного вектора. "Дорого" будет только по времени...
"нечист"
Многие в ХIIХ веке считали, что преобразование Фурье от "нечистого". Сам же автор свято верил в "Божественную гармонию"...

[Бахурин Сергей]

Ну во первых. Вы никогда не получите бпф на 2^27 точек через БПФ 2^16. В смысле конечно можно это сделать, но потребуется этап объединения 2^11 кусков по 2^16 степени точек БПФ. При объединении у вас будет еще 2^11 степени наборов поворотных коэффициентов. В результате объединения у вас будет большая большая куча комплексных умножений, каждое из которых будет выполнятся с округлением, и эти округления будут накапливаться и что получится в результате совсем не ясно. То что я описал и есть собственно идея БПФ заменить одно длинное многими короткими и потом объединить результат. Вопрос то как раз и стоит - как в результате объединения накапливаются ошибки округления умножений? Я лично не проверял и не знаю как они накапливаются (проверять честно не хочу ибо это очень долго).
Во вторых: Жан Батист Фурье жил на стыке XVIII и XIX веков, а века XIIX я честно не знаю.
В третьих: Если я что-то не понимаю или в чем-то неуверен (таких вещей очень много), то я предпочитаю прежде чем писать на форуме эти вещи понять или проверить, и уже потом писать пост, чего и вам желаю. Правда иногда я заблуждаюсь и тем не менее пишу неверные вещи, или вещи в которых я не совсем уверен, но это редко.

[Ivan Karamazov]

Моя мысль в том , что необязательно брать N= 2^27. Взять 2^16 и много раз прокрутить.

Вы заметили, что в самом начале темы я писал, что юзаю FFTW? Вообще говоря, Инструмент. Единственная незадача -- последовательность вычислений отпределяется числом отсчетов и режимом тюнинга. При этом единственное, что гарантируется -- число операций всегда O(n*log(n)), даже если n -- простое число. Т.е. за исключением вырожденных случаев, когда n не раскладывается по степеням 2,3,5 и (7 или 11), последовательность вычислений непредсказуема, и лучшее, что можно сделать -- взять как можно больше разных сигналов длины, превышающей реальные потребности. И все равно ничего не доказать. С обычным ДПФ (2**m) в этом смыле, кажется проще, но ненамного.
И это, заметьте не offtopic, а полезность, которую важно понимать всем, кто пользуется FFTW.

Сам же автор свято верил в "Божественную гармонию"...

"Чистая" математика, особенно времен Фурье -- провоцирует веру в высокое... А вот плавающая точка очень сильно развивает агностицизм -- это, конечно, offtopic.
----
А изложенная мною метода верификации меня практически устраивает. И понять ее не слишком сложно. А вот улучшить -- наверное можно.

[Santik]

Бахурин Сергей

Вы никогда не получите бпф на 2^27 точек через БПФ 2^16.

А это и не предполагалось!
Берем вещественный вектор длиной 2^16 + 666 (чисто из уважения к Инструменту) . Делаем FFT-> IFFT. Выходной вектор (комплексный) снова FFT-> IFFT и т.д.
После каждого этапа смотрим, как и предлагает уважаемый Ivan Karamazov :
(1) статистику разности действительной части комплексного сигнала и исходного сигнала;
(2) статистику подавления отрицательных частот в спектре комплексного сигнала (т.е. для этого еще одно ДПФ делается).

А вот плавающая точка очень сильно развивает агностицизм

Надежда всегда есть! А как же элементарная (для многих!) задача 0.5 л /3 - и без остатка!!!
-- это, конечно, offtopic.

[Бахурин Сергей]

и что же это даст? Если вы что в результате fft - ifft средний уровень мнимой части 1e-15 и что это значит? И как это соотнести с уровнем мнимой части при N = 2^27 ?

[Santik]

прежде чем писать на форуме эти вещи понять или проверить, и уже потом писать пост, чего и вам желаю

Ну вот, проверил:
do 10, i=1,n
x(i)=(0.,0.)
FX(i)=(0.,0.)
10 continue
M=100000000
n=8
x(2)=(1.,0.)
do 3 i=1,M
call FFTCF(n,X,FX)
call FFTCB(n,FX,X)
do 3 j=1,n
x(j)=x(j)/n
3 continue


M=100000000

0.000E+00 0.000E+00
0.100E+01 0.250E+00
0.000E+00 0.000E+00
0.313E-01 -0.245E+00
0.000E+00 0.000E+00
0.596E+00 0.240E+00
0.000E+00 0.000E+00
-0.312E-01 -0.245E+00

Compaq Visual FORTRAN 6.6.0
IMSL(R) Fortran 90 MP Library Version 4.01

[Бахурин Сергей]

Я не знаю что у вас за чудо фортран. Провел тот же эксперимент на своей длл-ке получил отличные от вашего фортрана результаты. Вот программа экперимента

Код: Выделить всё

#include < stdio.h >
#include < stdlib.h >
#include < windows.h >
#define _USE_MATH_DEFINES
#include < math.h >
#include "dsp.h"

HMODULE hDLL;



int main(){
	//загрузка dsp.dll
	hDLL = LoadDSP();
	if(hDLL)
		printf("Load DSP.DLL...............OK\n");
	else{
		printf("Load DSP.DLL...............ERROR\n");
		system("Pause");
		return 0;
	}
	
	int N = 8;
	//выделяю память под массивы
	double* sr = (double*)malloc(N*sizeof(double));
	double* si = (double*)malloc(N*sizeof(double));
	double* SR = (double*)malloc(N*sizeof(double));
	double* SI = (double*)malloc(N*sizeof(double));
	memset(sr,0,N*sizeof(double));
	memset(si,0,N*sizeof(double));
	sr[1] = 1.0;

	//создаю fft
	void* fft = NULL;
	fftCreate(N,&fft);
	
	int M = 100000000; // столько раз крутимся
	
	//цикл
	for(int i = 0; i<M; i++){
		if((i%(M/100)) == (M/100-1)) 
			printf("\r %d %%", (int)(100.0*(double)i/(double)M)+1); //вывожу % обработанного

		fftComplex(sr,si,SR,SI,N,fft); // комплексное fft
		ifftComplex(SR,SI,sr,si,N,fft);// комплексное ifft
	}
	printf("\n\n");

	double sum = 0.0; // cумма модуля мнимой части результата

	for(int i = 0; i<N; i++){
		sum+=fabs(si[i]); // накапливаю сумму
		printf("%+.12e\t%+.12e\n", sr[i], si[i]); // вывожу результат
	}
	printf("\n\n");
	printf("mean imag = %+.12e\n", sum/(double)N); // вывожу среднее мнимой части
	printf("\n\n");

	free(sr);
	free(si);
	free(SR);
	free(SI);
	fftClear(&fft);

	FreeLibrary(hDLL);
	system("Pause");
	return 0;
} 

dsp.dll можно скачать с сайта, там же dsp.h и dsp.cpp.
результаты:
Load DSP.DLL...............OK
100 %

+0.000000000000e+000 -0.000000000000e+000
+1.000000000000e+000 -0.000000000000e+000
+0.000000000000e+000 -0.000000000000e+000
-7.654042494671e-018 -7.654042494671e-018
+0.000000000000e+000 -0.000000000000e+000
+0.000000000000e+000 -0.000000000000e+000
+0.000000000000e+000 -0.000000000000e+000
+7.654042494671e-018 +7.654042494671e-018


mean imag = +1.913510623668e-018


Для продолжения нажмите любую клавишу . . .

ну как бы ничего не разваливается
пробовал N = 1024 M = 100000 получил
mean imag = +9.103270703708e-014
Тоже не разваливается

проверял при M = 1 N = 2^24 получил
mean imag = +1.102185690381e-020
и тоже не разваливается. Для больших длин проверить не смог ибо разрядности не хватает (система 32 бита) вываливается по эксепшену. Короче различаются наши алгоритмы бпф.