реклама на сайте
подробности

 
 
3 страниц V  < 1 2 3 >  
Reply to this topicStart new topic
> Быстрая свертка! Как?, помогите, пожалуйста, с алгоритмом
coolibin
сообщение Oct 15 2007, 13:36
Сообщение #16


Местный
***

Группа: Участник
Сообщений: 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 ищет ближайшую степень двойки которая больше чем заданое число


--------------------
Нет повести печальнее на свете, чем повесть о хреновом интернете.
Go to the top of the page
 
+Quote Post
shasik
сообщение Oct 15 2007, 13:58
Сообщение #17


Местный
***

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



В Matlab'e приведенный выше пример полностью прокатывает, поэтому еще раз обращаю внимание на вычисление прямого и обратного БПФ.
А еще смущает строка:
Код
for(i = 0; i != uSizeA; i++)


Если будете делать фильтрацию в лоб и с помощью БПФ, а затем посмотрите графики выходных отфильтрованных сигналов, то у Вас может наблюдаться небольшой сдвиг (в зависимости от того, как Вы поступили с начальными условиями - считали, что до этого сигнал был равен 0 либо подождали, пока не прибудет необходимое количество отсчетов). И если для сравнения результатов Вы вычисляли СКО сигналов - результаты будут не совпадать
Go to the top of the page
 
+Quote Post
coolibin
сообщение Oct 16 2007, 06:36
Сообщение #18


Местный
***

Группа: Участник
Сообщений: 214
Регистрация: 19-07-07
Пользователь №: 29 228



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


--------------------
Нет повести печальнее на свете, чем повесть о хреновом интернете.
Go to the top of the page
 
+Quote Post
shasik
сообщение Oct 18 2007, 08:05
Сообщение #19


Местный
***

Группа: Свой
Сообщений: 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 работает правильно!

Ну, тогда ладно, уговорили! Сдаемся! Говорите отгадку - что неправильно?

ЗЫ: расшифруйте понятие "неправильно" (архивчик с тестовым сигналом; с тем, что должно было получиться; с тем, что получилось на сам деле и с отсчетами фильтра очень помог бы).
Go to the top of the page
 
+Quote Post
coolibin
сообщение Oct 18 2007, 08:21
Сообщение #20


Местный
***

Группа: Участник
Сообщений: 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;
         }
     }
}


--------------------
Нет повести печальнее на свете, чем повесть о хреновом интернете.
Go to the top of the page
 
+Quote Post
shasik
сообщение Oct 22 2007, 10:00
Сообщение #21


Местный
***

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



Цитата(coolibin @ Oct 16 2007, 09:36) *
Алгоритм описаный здесь http://alglib.sources.ru/fft/fastconvolution.php


Оно Вам надо?
Я гораздо больше уважаю Алгоритмы для Программеров. Там можно скачать книгу (на 3 МБ в архиве) и сишные исходники (около 1 МБ). Почитав главу относящуюся к преобразованию Фурье и его применение, Вы многое поймете. Например, поймете, что зря использовали для своего алгоритма название fft, т.к. он далеко не быстрый. Может быть поэтому у Вас долго считает.

А на счет разделения на + и - части - реально с этим никто не заморачивается, все делается гороздо проще. Сходите по ссылке, рекомендую
Go to the top of the page
 
+Quote Post
coolibin
сообщение Oct 23 2007, 07:03
Сообщение #22


Местный
***

Группа: Участник
Сообщений: 214
Регистрация: 19-07-07
Пользователь №: 29 228



А никто не скажет, когда применять быструю свертку, а когда обычную? или всегда выгодно быструю?


--------------------
Нет повести печальнее на свете, чем повесть о хреновом интернете.
Go to the top of the page
 
+Quote Post
shasik
сообщение Oct 23 2007, 07:26
Сообщение #23


Местный
***

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



Цитата(coolibin @ Oct 23 2007, 10:03) *
А никто не скажет, когда применять быструю свертку, а когда обычную? или всегда выгодно быструю?

Все зависит: от длины используемого фильтра (скользящее среднее чаще всего делается напрямую без использования быстрых алгоритмов), от используемого железа (многие процессоры ЦОС поддерживают выполнение свертки на аппаратном уровне причем параллельно с основным железом), от требования к быстродействию программы, от наличия времени на написание программы (свертка в лоб всегда легче пишется), от того доступен ли весь сигнал сразу или накопление его отсчетов идет в real-time да и расположение свезд на небе тоже наверное влияет каким-нибудь образом.

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

ЗЫ: Еще раз говорю: скачайте то, что я рекомендовал! Там есть, что нужно.
Go to the top of the page
 
+Quote Post
coolibin
сообщение Oct 23 2007, 08:00
Сообщение #24


Местный
***

Группа: Участник
Сообщений: 214
Регистрация: 19-07-07
Пользователь №: 29 228



Цитата(shasik @ Oct 23 2007, 10:26) *
ЗЫ: Еще раз говорю: скачайте то, что я рекомендовал! Там есть, что нужно.


А что именно? там много всякой всячины!

Книга fxtbook.pdf что ли?

Сообщение отредактировал coolibin - Oct 23 2007, 08:13


--------------------
Нет повести печальнее на свете, чем повесть о хреновом интернете.
Go to the top of the page
 
+Quote Post
fontp
сообщение Oct 23 2007, 08:19
Сообщение #25


Эксперт
*****

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



Цитата(coolibin @ Oct 23 2007, 12:00) *
А что именно? там много всякой всячины!

Книга fxtbook.pdf что ли?


Всё что описано в книге - реализовано в fxt package на C
http://www.jjj.de/fxt/fxtpage.html#fxt
Go to the top of the page
 
+Quote Post
shasik
сообщение Oct 23 2007, 11:43
Сообщение #26


Местный
***

Группа: Свой
Сообщений: 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

Ксати, специально посмотрел соответствующий раздел - Ваша мечта сбылась, там тоже есть разделение на +/- часть.
ЗЫ: надеюсь переведете с английского уже без посторонней помощи
Go to the top of the page
 
+Quote Post
Ole2
сообщение Dec 3 2010, 15:54
Сообщение #27





Группа: Новичок
Сообщений: 5
Регистрация: 2-12-10
Пользователь №: 61 351



Подскажите пожалуйста, откуда берется выигрыш при проведении свертки или корреляции в частотной области? Количество умножений и суммирований одинаковое в частотной и временной областях. Помимо этого для частотной области добавляется прямое и обратное преобразование Фурье. Получается проигрыш. Гмм...
Go to the top of the page
 
+Quote Post
fontp
сообщение Dec 3 2010, 20:43
Сообщение #28


Эксперт
*****

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



QUOTE (Ole2 @ Dec 3 2010, 18:54) *
Подскажите пожалуйста, откуда берется выигрыш при проведении свертки или корреляции в частотной области? Количество умножений и суммирований одинаковое в частотной и временной областях. Помимо этого для частотной области добавляется прямое и обратное преобразование Фурье. Получается проигрыш. Гмм...


Гмм...Где ж одинаковое? В частотной области свертка (N*K операций) превращается просто в произведение спектров (N).
Go to the top of the page
 
+Quote Post
Ole2
сообщение Dec 4 2010, 18:54
Сообщение #29





Группа: Новичок
Сообщений: 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. Чудовищный проигрыш…
Go to the top of the page
 
+Quote Post
fontp
сообщение Dec 4 2010, 19:41
Сообщение #30


Эксперт
*****

Группа: Свой
Сообщений: 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 может быть предпочтительней только если число сдвигов К мало
Go to the top of the page
 
+Quote Post

3 страниц V  < 1 2 3 >
Reply to this topicStart new topic
1 чел. читают эту тему (гостей: 1, скрытых пользователей: 0)
Пользователей: 0

 


RSS Текстовая версия Сейчас: 6th July 2025 - 11:15
Рейтинг@Mail.ru


Страница сгенерированна за 0.0153 секунд с 7
ELECTRONIX ©2004-2016