|
|
  |
Быстрая свертка! Как?, помогите, пожалуйста, с алгоритмом |
|
|
|
Oct 15 2007, 13:36
|
Местный
  
Группа: Участник
Сообщений: 214
Регистрация: 19-07-07
Пользователь №: 29 228

|
Цитата(shasik @ Oct 15 2007, 13:08)  Почему Код uint uConvSize = get_length(uSizeA + uSizeB - 1); а не просто Код uint uConvSize = uSizeA + uSizeB - 1; Функция get_length ищет ближайшую степень двойки которая больше чем заданое число
--------------------
Нет повести печальнее на свете, чем повесть о хреновом интернете.
|
|
|
|
|
Oct 15 2007, 13:58
|

Местный
  
Группа: Свой
Сообщений: 319
Регистрация: 3-09-05
Из: Беларусь, Новополоцк
Пользователь №: 8 188

|
В Matlab'e приведенный выше пример полностью прокатывает, поэтому еще раз обращаю внимание на вычисление прямого и обратного БПФ. А еще смущает строка: Код for(i = 0; i != uSizeA; i++) Если будете делать фильтрацию в лоб и с помощью БПФ, а затем посмотрите графики выходных отфильтрованных сигналов, то у Вас может наблюдаться небольшой сдвиг (в зависимости от того, как Вы поступили с начальными условиями - считали, что до этого сигнал был равен 0 либо подождали, пока не прибудет необходимое количество отсчетов). И если для сравнения результатов Вы вычисляли СКО сигналов - результаты будут не совпадать
|
|
|
|
|
Oct 16 2007, 06:36
|
Местный
  
Группа: Участник
Сообщений: 214
Регистрация: 19-07-07
Пользователь №: 29 228

|
Алгоритм описаный здесь http://alglib.sources.ru/fft/fastconvolution.php, т .е. отклик делится на две части полож. и отриц., работаем прекрасно и на моем БПФ. Значит не в БПФ дело
--------------------
Нет повести печальнее на свете, чем повесть о хреновом интернете.
|
|
|
|
|
Oct 18 2007, 08:05
|

Местный
  
Группа: Свой
Сообщений: 319
Регистрация: 3-09-05
Из: Беларусь, Новополоцк
Пользователь №: 8 188

|
Цитата(coolibin @ Oct 12 2007, 14:39)  Что в этой свертке неправильно? Код void fast_conv(const double *pSigA, unsigned int uSizeA, const double *pSigB,unsigned int uSizeB, double *pSigR){ complex<double> *pSigACx, *pSigBCx;
.......................................................................
delete [] pSigACx; delete [] pSigBCx; } fft работает правильно! Ну, тогда ладно, уговорили! Сдаемся! Говорите отгадку - что неправильно? ЗЫ: расшифруйте понятие "неправильно" (архивчик с тестовым сигналом; с тем, что должно было получиться; с тем, что получилось на сам деле и с отсчетами фильтра очень помог бы).
|
|
|
|
|
Oct 18 2007, 08:21
|
Местный
  
Группа: Участник
Сообщений: 214
Регистрация: 19-07-07
Пользователь №: 29 228

|
Цитата(shasik @ Oct 18 2007, 11:05)  Ну, тогда ладно, уговорили! Сдаемся! Говорите отгадку - что неправильно?
ЗЫ: расшифруйте понятие "неправильно" (архивчик с тестовым сигналом; с тем, что должно было получиться; с тем, что получилось на сам деле и с отсчетами фильтра очень помог бы). Я не знаю что неправильно. Дело в том, что метод разделения отклика на положительный и отрицательный работает и на моем БПФ и мне этого достаточно. Вот рабочая свертка: Код //быстрая свертка void COOLibin_math::fastconvolution(std::vector<double>& Signal, const std::vector<double>& Responce, uint uPositiveLen) { std::vector<std::complex<double> > SignalCx, ResponceCx; uint uConvLen, uSignalLen, uNegativeLen; size_t i;
uSignalLen = Signal.size(); uNegativeLen = Responce.size() - uPositiveLen - 1;
uConvLen = uSignalLen; if(uNegativeLen > uPositiveLen ){ uConvLen += uNegativeLen; } else{ uConvLen += uPositiveLen; }
if( uNegativeLen + uPositiveLen + 1 > uConvLen ){ uConvLen = uNegativeLen + uPositiveLen + 1; }
i = 1; while(i < uConvLen) { i = i*2; } uConvLen = i;
SignalCx.resize(uConvLen); ResponceCx.resize(uConvLen);
for(i = 0; i != uSignalLen; i++){ SignalCx[i] = std::complex<double>(Signal[i], 0.0); }
for(i = uSignalLen; i != uConvLen; i++){ SignalCx[i] = std::complex<double>(0.0, 0.0); }
for(i = 0; i != uConvLen; i++){ ResponceCx[i] = std::complex<double>(0.0, 0.0); }
for(i = 0; i <= uPositiveLen; i++){ ResponceCx[i] = std::complex<double>(Responce[i + uNegativeLen], 0.0); }
for(i = 1; i <= uNegativeLen; i++){ ResponceCx[uConvLen - i] = std::complex<double>(Responce[uNegativeLen - i], 0.0); }
fft(SignalCx, false); fft(ResponceCx, false);
for(i = 0; i != uConvLen; i++){ SignalCx[i] *= ResponceCx[i]; }
fft(SignalCx, true);
for(i = 0; i != uSignalLen; i++){ Signal[i] = SignalCx[i].real(); } } А это сам БПФ: Код void COOLibin_math::rec_fft(std::vector<std::complex<double> >& Signal, bool bInv){
uint uSignalLen = Signal.size(); if(uSignalLen == 1){ return; }
const double PI = 3.1415926535897932; std::vector<std::complex<double> > Odd(uSignalLen/2), Even(uSignalLen/2);
for(size_t i = 0; i != uSignalLen/2; i++){ Odd[i] = Signal[2*i]; Even[i] = Signal[2*i + 1]; }
rec_fft(Odd, bInv); rec_fft(Even, bInv);
double fSign; if(!bInv){ fSign = 1.0; } else{ fSign = -1.0; }
std::complex<double> W(1.0, 0.0), WStep(cos(2.0*PI/uSignalLen), sin(fSign*2.0*PI/uSignalLen)), Temp;
for(size_t k = 0; k != uSignalLen/2; k++){ Temp = W*Even[k]; Signal[k] = Odd[k] + Temp; Signal[k + uSignalLen/2] = Odd[k] - Temp; W *= WStep; } }
void COOLibin_math::fft(std::vector<std::complex<double> >& Signal, bool bInv){
rec_fft(Signal, bInv);
if(bInv){ uint uSignalLen = Signal.size(); for(std::vector<std::complex<double> >::iterator it = Signal.begin(); it != Signal.end(); it++){ (*it) /= uSignalLen; } } }
--------------------
Нет повести печальнее на свете, чем повесть о хреновом интернете.
|
|
|
|
|
Oct 22 2007, 10:00
|

Местный
  
Группа: Свой
Сообщений: 319
Регистрация: 3-09-05
Из: Беларусь, Новополоцк
Пользователь №: 8 188

|
Цитата(coolibin @ Oct 16 2007, 09:36)  Оно Вам надо? Я гораздо больше уважаю Алгоритмы для Программеров. Там можно скачать книгу (на 3 МБ в архиве) и сишные исходники (около 1 МБ). Почитав главу относящуюся к преобразованию Фурье и его применение, Вы многое поймете. Например, поймете, что зря использовали для своего алгоритма название fft, т.к. он далеко не быстрый. Может быть поэтому у Вас долго считает. А на счет разделения на + и - части - реально с этим никто не заморачивается, все делается гороздо проще. Сходите по ссылке, рекомендую
|
|
|
|
|
Oct 23 2007, 07:26
|

Местный
  
Группа: Свой
Сообщений: 319
Регистрация: 3-09-05
Из: Беларусь, Новополоцк
Пользователь №: 8 188

|
Цитата(coolibin @ Oct 23 2007, 10:03)  А никто не скажет, когда применять быструю свертку, а когда обычную? или всегда выгодно быструю? Все зависит: от длины используемого фильтра (скользящее среднее чаще всего делается напрямую без использования быстрых алгоритмов), от используемого железа (многие процессоры ЦОС поддерживают выполнение свертки на аппаратном уровне причем параллельно с основным железом), от требования к быстродействию программы, от наличия времени на написание программы (свертка в лоб всегда легче пишется), от того доступен ли весь сигнал сразу или накопление его отсчетов идет в real-time да и расположение свезд на небе тоже наверное влияет каким-нибудь образом. А если серьезно то в первом приближении оценивается по количеству операций - на малых длинах быстрые алгоритмы дают совсем незначительный выйгрыш по быстродействию, а иногда даже медленнее решения в лоб, а гемороя от них много. ЗЫ: Еще раз говорю: скачайте то, что я рекомендовал! Там есть, что нужно.
|
|
|
|
|
Oct 23 2007, 08:00
|
Местный
  
Группа: Участник
Сообщений: 214
Регистрация: 19-07-07
Пользователь №: 29 228

|
Цитата(shasik @ Oct 23 2007, 10:26)  ЗЫ: Еще раз говорю: скачайте то, что я рекомендовал! Там есть, что нужно. А что именно? там много всякой всячины! Книга fxtbook.pdf что ли?
Сообщение отредактировал coolibin - Oct 23 2007, 08:13
--------------------
Нет повести печальнее на свете, чем повесть о хреновом интернете.
|
|
|
|
|
Oct 23 2007, 11:43
|

Местный
  
Группа: Свой
Сообщений: 319
Регистрация: 3-09-05
Из: Беларусь, Новополоцк
Пользователь №: 8 188

|
Цитата(coolibin @ Oct 23 2007, 11:00)  А что именно? там много всякой всячины! Книга fxtbook.pdf что ли? Ну, Вы и ленивый! Книга в формате pdf в архиве - http://www.jjj.de/fxt/fxtbook.pdf.gzКнига в формате dvi в архиве - http://www.jjj.de/fxt/fxtbook.dvi.gz - самый небольшой размер Исходники (также в архиве) - http://www.jjj.de/fxt/fxt-2007.10.14.tgzКсати, специально посмотрел соответствующий раздел - Ваша мечта сбылась, там тоже есть разделение на +/- часть. ЗЫ: надеюсь переведете с английского уже без посторонней помощи
|
|
|
|
|
Dec 3 2010, 15:54
|
Группа: Новичок
Сообщений: 5
Регистрация: 2-12-10
Пользователь №: 61 351

|
Подскажите пожалуйста, откуда берется выигрыш при проведении свертки или корреляции в частотной области? Количество умножений и суммирований одинаковое в частотной и временной областях. Помимо этого для частотной области добавляется прямое и обратное преобразование Фурье. Получается проигрыш. Гмм...
|
|
|
|
|
Dec 4 2010, 18:54
|
Группа: Новичок
Сообщений: 5
Регистрация: 2-12-10
Пользователь №: 61 351

|
Ну ладно… Пусть мы имеем отсчеты квадратурного (комплексного) сигнала, поступающие непрерывно на коррелятор, имеющий линию задержки длинной k с комплексными эталонными коэффициентами. Один такт вычислений в этом случае даст 4 операции умножения и 3 суммирования на 1 ячейку линии задержки (комплексное умножение). Так как этих ячеек k, то общее количество операций k*(4умн+3сум). Кроме того, необходио суммирование результата умножения по всем ячейкам. Итого количество операций Ot для корреляция во времени: Ot= k*(4умн+3сум)+kсум; Для частотной области только перемножаем спектры сигналов и получаем количество опеаций k*(4умн+3сум). Суммировать здесь действительно не надо. Но для выполнения БПФ необходимо k*log2(k) операций. Итого для вычисления свертки в частотной области необходимо Of операций: Of= k*(4умн+3сум)+2*k*log2(k). Так при k=2048 для вычисления свертки во времени необходимо 16384 операций умножения и суммирования, а для частоной 59392. Чудовищный проигрыш…
|
|
|
|
|
Dec 4 2010, 19:41
|

Эксперт
    
Группа: Свой
Сообщений: 1 467
Регистрация: 25-06-04
Пользователь №: 183

|
QUOTE (Ole2 @ Dec 4 2010, 21:54)  Ну ладно… Пусть мы имеем отсчеты квадратурного (комплексного) сигнала, поступающие непрерывно на коррелятор, имеющий линию задержки длинной k с комплексными эталонными коэффициентами. Один такт вычислений в этом случае даст 4 операции умножения и 3 суммирования на 1 ячейку линии задержки (комплексное умножение). Так как этих ячеек k, то общее количество операций k*(4умн+3сум). Кроме того, необходио суммирование результата умножения по всем ячейкам. Итого количество операций Ot для корреляция во времени: Ot= k*(4умн+3сум)+kсум; Для частотной области только перемножаем спектры сигналов и получаем количество опеаций k*(4умн+3сум). Суммировать здесь действительно не надо. Но для выполнения БПФ необходимо k*log2(k) операций. Итого для вычисления свертки в частотной области необходимо Of операций: Of= k*(4умн+3сум)+2*k*log2(k). Так при k=2048 для вычисления свертки во времени необходимо 16384 операций умножения и суммирования, а для частоной 59392. Чудовищный проигрыш… Ерунда... Свертка это сумма x(i)*x(i+k) длинною N, для всех к от 0 до K. Всего N*K операций, порядка N*N в ассимптотике, K<N 2 фурье это порядка 2*N*log(N) плюс N умножений спектров для всех K сразу Вот и сравните N*log(N) и N*N при больших N Прямая свертка при больших N может быть предпочтительней только если число сдвигов К мало
|
|
|
|
|
  |
1 чел. читают эту тему (гостей: 1, скрытых пользователей: 0)
Пользователей: 0
|
|
|