Устойчивость БПФ/ОБПФ на "длинных" сигналах
-
- Сообщения: 89
- Зарегистрирован: 28 окт 2010, 22:31
- Откуда: Москва
Re: Устойчивость БПФ/ОБПФ на "длинных" сигналах
Про последний тезис -- более-менне в курсе. Речь только о возможном "развале" длинного сигнала из-за ошибок округления при ДПФ.
А "насчет зависимости от характера сигнала" -- все оказалось очень просто (как и предполагалось) -- в тестовом сигнале присутствовало одно локальное ограничение ("clip") на 4-5 остсчетов. Это дало "выбоину" на одном отсчете (причем только у мнимой части) Гильберта. Т.е. не фатальную выбоину -- просто на несколько младших разрядов у разных реализаций БПФ. (Ну да, избегайте ограничений -- они имеют "богатый" спектр со всеми вытекющими). И, как уже упоминалось -- если ваш сигнал имеет четное число отсчетов (при необходимости его легко разумным образом "доопределить") -- необходимо обнулить в т.ч. и n/2-ю компоненту спектра.
А "насчет зависимости от характера сигнала" -- все оказалось очень просто (как и предполагалось) -- в тестовом сигнале присутствовало одно локальное ограничение ("clip") на 4-5 остсчетов. Это дало "выбоину" на одном отсчете (причем только у мнимой части) Гильберта. Т.е. не фатальную выбоину -- просто на несколько младших разрядов у разных реализаций БПФ. (Ну да, избегайте ограничений -- они имеют "богатый" спектр со всеми вытекющими). И, как уже упоминалось -- если ваш сигнал имеет четное число отсчетов (при необходимости его легко разумным образом "доопределить") -- необходимо обнулить в т.ч. и n/2-ю компоненту спектра.
Если ваши решения вам нравятся -- это хорошие решения. И наоборот.
- Santik
- Сообщения: 609
- Зарегистрирован: 28 дек 2010, 08:04
- Откуда: Мирный (Якутия)
- Контактная информация:
Re: Устойчивость БПФ/ОБПФ на "длинных" сигналах
А сколько пар преобразований (FFT->IFFT) -> (FFT->IFFT) ->...можно сделать без "развала" исходного сигнала длиной N?
Можно оценить Nmax=N*2^M
Можно оценить Nmax=N*2^M
-
- Сообщения: 89
- Зарегистрирован: 28 окт 2010, 22:31
- Откуда: Москва
Re: Устойчивость БПФ/ОБПФ на "длинных" сигналах
Вообще, такой эксперимент дорог и нечист, хотя, конечно, попробовать можно. На самом деле я сделал так -- ввел режим верификации, выводящий
(1) статистику разности действительной части комплексного сигнала и исходного сигнала;
(2) статистику подавления отрицательных частот в спектре комплексного сигнала (т.е. для этого еще одно ДПФ делается).
Результаты -- разница в действительных частях в пределах младшего (из исходных 16-ти) бит, подавление -- не хуже -170 дБ (обычно порядка 200 (!) дБ. Но это -- дорого, поскольку у меня всего 4 ГБ RAM и Win7, у которой после однократного перерасхода памяти любит возникать ситуация: мой процесс в данный момент потребляет всего 2,7 ГБ, а размер дискового кэша -- 3,5 ГБ! Естественно, все стоит, а кэш общается с диском. (Такую систему угробили )
(1) статистику разности действительной части комплексного сигнала и исходного сигнала;
(2) статистику подавления отрицательных частот в спектре комплексного сигнала (т.е. для этого еще одно ДПФ делается).
Результаты -- разница в действительных частях в пределах младшего (из исходных 16-ти) бит, подавление -- не хуже -170 дБ (обычно порядка 200 (!) дБ. Но это -- дорого, поскольку у меня всего 4 ГБ RAM и Win7, у которой после однократного перерасхода памяти любит возникать ситуация: мой процесс в данный момент потребляет всего 2,7 ГБ, а размер дискового кэша -- 3,5 ГБ! Естественно, все стоит, а кэш общается с диском. (Такую систему угробили )
Если ваши решения вам нравятся -- это хорошие решения. И наоборот.
- Santik
- Сообщения: 609
- Зарегистрирован: 28 дек 2010, 08:04
- Откуда: Мирный (Якутия)
- Контактная информация:
Re: Устойчивость БПФ/ОБПФ на "длинных" сигналах
Извините, опять "по диагонали" все прочитал. Но чтобы Вашу методику осознать мне потребуется очень много времени.Вообще, такой эксперимент дорог и нечист
"дорог"
Моя мысль в том , что необязательно брать N= 2^27. Взять 2^16 и много раз прокрутить. Памяти много не надо на это. И следить за мнимой частью выходного вектора. "Дорого" будет только по времени...
"нечист"
Многие в ХIIХ веке считали, что преобразование Фурье от "нечистого". Сам же автор свято верил в "Божественную гармонию"...
- Бахурин Сергей
- Администратор
- Сообщения: 1116
- Зарегистрирован: 05 окт 2010, 19:55
- Контактная информация:
Re: Устойчивость БПФ/ОБПФ на "длинных" сигналах
Ну во первых. Вы никогда не получите бпф на 2^27 точек через БПФ 2^16. В смысле конечно можно это сделать, но потребуется этап объединения 2^11 кусков по 2^16 степени точек БПФ. При объединении у вас будет еще 2^11 степени наборов поворотных коэффициентов. В результате объединения у вас будет большая большая куча комплексных умножений, каждое из которых будет выполнятся с округлением, и эти округления будут накапливаться и что получится в результате совсем не ясно. То что я описал и есть собственно идея БПФ заменить одно длинное многими короткими и потом объединить результат. Вопрос то как раз и стоит - как в результате объединения накапливаются ошибки округления умножений? Я лично не проверял и не знаю как они накапливаются (проверять честно не хочу ибо это очень долго).
Во вторых: Жан Батист Фурье жил на стыке XVIII и XIX веков, а века XIIX я честно не знаю.
В третьих: Если я что-то не понимаю или в чем-то неуверен (таких вещей очень много), то я предпочитаю прежде чем писать на форуме эти вещи понять или проверить, и уже потом писать пост, чего и вам желаю. Правда иногда я заблуждаюсь и тем не менее пишу неверные вещи, или вещи в которых я не совсем уверен, но это редко.
Во вторых: Жан Батист Фурье жил на стыке XVIII и XIX веков, а века XIIX я честно не знаю.
В третьих: Если я что-то не понимаю или в чем-то неуверен (таких вещей очень много), то я предпочитаю прежде чем писать на форуме эти вещи понять или проверить, и уже потом писать пост, чего и вам желаю. Правда иногда я заблуждаюсь и тем не менее пишу неверные вещи, или вещи в которых я не совсем уверен, но это редко.
-
- Сообщения: 89
- Зарегистрирован: 28 окт 2010, 22:31
- Откуда: Москва
Re: Устойчивость БПФ/ОБПФ на "длинных" сигналах
Вы заметили, что в самом начале темы я писал, что юзаю FFTW? Вообще говоря, Инструмент. Единственная незадача -- последовательность вычислений отпределяется числом отсчетов и режимом тюнинга. При этом единственное, что гарантируется -- число операций всегда O(n*log(n)), даже если n -- простое число. Т.е. за исключением вырожденных случаев, когда n не раскладывается по степеням 2,3,5 и (7 или 11), последовательность вычислений непредсказуема, и лучшее, что можно сделать -- взять как можно больше разных сигналов длины, превышающей реальные потребности. И все равно ничего не доказать. С обычным ДПФ (2**m) в этом смыле, кажется проще, но ненамного.Моя мысль в том , что необязательно брать N= 2^27. Взять 2^16 и много раз прокрутить.
И это, заметьте не offtopic, а полезность, которую важно понимать всем, кто пользуется FFTW.
"Чистая" математика, особенно времен Фурье -- провоцирует веру в высокое... А вот плавающая точка очень сильно развивает агностицизм -- это, конечно, offtopic.Сам же автор свято верил в "Божественную гармонию"...
----
А изложенная мною метода верификации меня практически устраивает. И понять ее не слишком сложно. А вот улучшить -- наверное можно.
Если ваши решения вам нравятся -- это хорошие решения. И наоборот.
- Santik
- Сообщения: 609
- Зарегистрирован: 28 дек 2010, 08:04
- Откуда: Мирный (Якутия)
- Контактная информация:
Re: Устойчивость БПФ/ОБПФ на "длинных" сигналах
Бахурин Сергей
Берем вещественный вектор длиной 2^16 + 666 (чисто из уважения к Инструменту) . Делаем FFT-> IFFT. Выходной вектор (комплексный) снова FFT-> IFFT и т.д.
После каждого этапа смотрим, как и предлагает уважаемый Ivan Karamazov :
(1) статистику разности действительной части комплексного сигнала и исходного сигнала;
(2) статистику подавления отрицательных частот в спектре комплексного сигнала (т.е. для этого еще одно ДПФ делается).
-- это, конечно, offtopic.
А это и не предполагалось!Вы никогда не получите бпф на 2^27 точек через БПФ 2^16.
Берем вещественный вектор длиной 2^16 + 666 (чисто из уважения к Инструменту) . Делаем FFT-> IFFT. Выходной вектор (комплексный) снова FFT-> IFFT и т.д.
После каждого этапа смотрим, как и предлагает уважаемый Ivan Karamazov :
(1) статистику разности действительной части комплексного сигнала и исходного сигнала;
(2) статистику подавления отрицательных частот в спектре комплексного сигнала (т.е. для этого еще одно ДПФ делается).
Надежда всегда есть! А как же элементарная (для многих!) задача 0.5 л /3 - и без остатка!!!А вот плавающая точка очень сильно развивает агностицизм
-- это, конечно, offtopic.
- Бахурин Сергей
- Администратор
- Сообщения: 1116
- Зарегистрирован: 05 окт 2010, 19:55
- Контактная информация:
Re: Устойчивость БПФ/ОБПФ на "длинных" сигналах
и что же это даст? Если вы что в результате fft - ifft средний уровень мнимой части 1e-15 и что это значит? И как это соотнести с уровнем мнимой части при N = 2^27 ?
- Santik
- Сообщения: 609
- Зарегистрирован: 28 дек 2010, 08:04
- Откуда: Мирный (Якутия)
- Контактная информация:
Re: Устойчивость БПФ/ОБПФ на "длинных" сигналах
Ну вот, проверил:прежде чем писать на форуме эти вещи понять или проверить, и уже потом писать пост, чего и вам желаю
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
- Бахурин Сергей
- Администратор
- Сообщения: 1116
- Зарегистрирован: 05 окт 2010, 19:55
- Контактная информация:
Re: Устойчивость БПФ/ОБПФ на "длинных" сигналах
Я не знаю что у вас за чудо фортран. Провел тот же эксперимент на своей длл-ке получил отличные от вашего фортрана результаты. Вот программа экперимента
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 бита) вываливается по эксепшену. Короче различаются наши алгоритмы бпф.
Код: Выделить всё
#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;
}
результаты:
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 бита) вываливается по эксепшену. Короче различаются наши алгоритмы бпф.