Форум программистов, компьютерный форум, киберфорум
Наши страницы
Assembler: математика, вычисления
Войти
Регистрация
Восстановить пароль
 
 
Рейтинг 4.94/34: Рейтинг темы: голосов - 34, средняя оценка - 4.94
vital792
2002 / 1274 / 60
Регистрация: 05.06.2010
Сообщений: 2,213
1

Преобразование Фурье на ассемблере

14.07.2010, 10:44. Просмотров 6275. Ответов 43
Метки нет (Все метки)

Ну так и здравствуйте! Имеется код на с (функция выполняет быстрое преобразование фурье):
C
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
void Fft2(double *in_Buffer, double *out_Buffer, long fftFrameSize, long sign)
{
    double wr, wi, arg, *p1, *p2, temp;
    double tr, ti, ur, ui, *p1r, *p1i, *p2r, *p2i;
    long i, bitm, j, le, le2, k;
 
    memcpy(out_Buffer, in_Buffer, fftFrameSize*2*sizeof(double));
 
    for (i = 2; i < 2*fftFrameSize-2; i += 2)
    {
        for (bitm = 2, j = 0; bitm < 2*fftFrameSize; bitm <<= 1)
        {
            if (i & bitm) j++;
            j <<= 1;
        }
        if (i < j)
        {
            p1 = out_Buffer+i; p2 = out_Buffer+j;
            temp = *p1; *(p1++) = *p2;
            *(p2++) = temp; temp = *p1;
            *p1 = *p2; *p2 = temp;
        }
    }
    long lim = (long)(log((double)fftFrameSize)/log(2.)+.5);
    for (k = 0, le = 2; k < lim; k++)
    {
        le <<= 1;
        le2 = le>>1;
        ur = 1.0;
        ui = 0.0;
        arg = pi / (le2>>1);
        wr = cos(arg);
        wi = sign*sin(arg);
        for (j = 0; j < le2; j += 2)
        {
            p1r = out_Buffer+j; p1i = p1r+1;
            p2r = p1r+le2; p2i = p2r+1;
            for (i = j; i < 2*fftFrameSize; i += le)
            {
                tr = *p2r * ur - *p2i * ui;
                ti = *p2r * ui + *p2i * ur;
                *p2r = *p1r - tr; *p2i = *p1i - ti;
                *p1r += tr; *p1i += ti;
                p1r += le; p1i += le;
                p2r += le; p2i += le;
            }
            tr = ur*wr - ui*wi;
            ui = ur*wi + ui*wr;
            ur = tr;
        }
    }
}
Я переписал его на ассемблер:
Assembler
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
.686p
.XMM
 
.model c,flat
 
; parametrs
in_data     equ dword ptr [ebp+8]
out_data    equ dword ptr [ebp+12]
size_data   equ dword ptr [ebp+16]
 
; local
dtemp   equ dword ptr [ebp-4]   ;dd 0
ftemp   equ qword ptr [ebp-12]  ;dq 0.0
wr  equ qword ptr [ebp-20]  ;dq 0.0
wi  equ qword ptr [ebp-28]  ;dq 0.0
 
 .code
 
public FFT_SSE2
 
; void FFT_SSE2( double *in_data, double *out_data, int fftFrameSize)
FFT_SSE2 proc
push    ebp
mov ebp, esp
;add    ebp, 8
sub esp, 28
 
; memcpy(out_Buffer, in_Buffer, fftFrameSize*2*sizeof(double));
 
mov ecx, size_data
shl ecx, 2              ; * 2 * sizeof(double)
mov esi, in_data
mov edi, out_data
rep movsd
 
mov ecx, 2              ; ecx - i
_l0:
 
mov edx, 2              ; edx - bitm
xor ebx, ebx            ; ebx - j
_l1:
 
mov eax, ecx
and eax, edx
jz  _l3
inc ebx
_l3:
shl ebx, 1
 
shl edx, 1
mov eax, size_data
shl eax, 1
cmp edx, eax
jb  _l1
 
cmp ecx, ebx
jnb _l4
mov eax, out_data
movapd  xmm0, qword ptr [eax+ecx*8] ; p1 = out_Buffer+i
movapd  xmm1, qword ptr [eax+ebx*8] ; p2 = out_Buffer+j
movapd  qword ptr [eax+ecx*8], xmm1
movapd  qword ptr [eax+ebx*8], xmm0
_l4:
 
add ecx, 2
mov eax, size_data
shl eax, 1
sub eax, 2
cmp ecx, eax
jb  _l0
 
fld1
fild    size_data
fyl2x
frndint
fistp   dtemp               ; (long)log2(size)
mov edx, dtemp
 
xor ecx, ecx                ; ecx - k
mov ebx, 2              ; ebx - le
main_loop:
 
push    edx
mov eax, ebx
shl ebx, 1
fldpi
shr eax, 1
mov dtemp, eax
fidiv   dtemp               ; st(0): arg = pi / (le2>>1)
fsincos                 ; st(0)=cos(arg) st(1)=sin(arg)
fstp    wr
fchs
fstp    wi
xorpd   xmm2, xmm2          ; ur|ui
fld1
fstp    ftemp
movlps  xmm2, qword ptr [ftemp]
 
xor esi, esi            ; esi - j
_loop1:
 
mov eax, out_data
mov edx, esi
shl edx, 3
add eax, edx            ; eax = out_data + j
 
mov edi, esi
_loop2:
 
movapd  xmm3, xmm2  ; push xmm2
 
movapd  xmm0, xmm2
movddup xmm1, qword ptr [eax+ebx*4+8]
movddup xmm2, qword ptr [eax+ebx*4]
mulpd   xmm2, xmm0
shufpd  xmm0, xmm0, 1
mulpd   xmm1, xmm0
addsubpd xmm2, xmm1
 
movapd  xmm0, xmmword ptr [eax]     ; *p1
subpd   xmm0, xmm2
movapd  xmmword ptr [eax+ebx*4], xmm0
 
movapd  xmm0, xmmword ptr [eax]
addpd   xmm0, xmm2
movapd  xmmword ptr [eax], xmm0
 
movapd  xmm2, xmm3  ; pop xmm2
 
push    ebx
shl ebx, 3
add eax, ebx
pop ebx
 
add edi, ebx
mov edx, size_data
shl edx, 1
cmp edi, edx
jb  _loop2
 
movapd  xmm0, xmm2
movddup xmm1, qword ptr [wi]
movddup xmm2, qword ptr [wr]
mulpd   xmm2, xmm0
shufpd  xmm0, xmm0, 1
mulpd   xmm1, xmm0
addsubpd xmm2, xmm1
 
add esi, 2
mov eax, ebx
shr eax, 1
cmp esi, eax
jb  _loop1
 
inc ecx
pop edx
cmp ecx, edx
jb  main_loop
 
mov esp, ebp
pop ebp
retn
FFT_SSE2 endp
 
end
Все работает как надо. Но я ожидал что использование sse3 увеличит производительность и функция будет работать хотя бы вдвое быстрее аналога с сопроцессором. В результате моя функция работает быстрее сишной тока на 15%. Посмотрите пожалуйста, можно ли как то оптимизировать - почему так медленно?
ЗЫ: память выровнена по адресам кратным 16
0
Similar
Эксперт
41792 / 34177 / 6122
Регистрация: 12.04.2006
Сообщений: 57,940
14.07.2010, 10:44
Ответы с готовыми решениями:

Найти коэффиценты разложения в ряд Фурье, используя быстрое преобразование Фурье (БПФ)
Прошу помочь мне в нелеггкой задачке нужно для заданной на периоде 2∏ функции...

Быстрое преобразование Фурье и ошибка "Неявное преобразование типа"
А подскажите еще по одной прблемке: Есть программа реализующая БПФ (ну должна...

Преобразование Фурье
Добрый день! Не могли бы вы подбросить рабочий код для преобразования Фурье...

Преобразование Фурье
К данному сигналу добавить шум y(t)=y(t)+rnd(1) Построить график с шумом И...

Преобразование Фурье
По форуму искал, подобные вопросы нашел, а вот ответов на них нету. Имеем...

43
D@rkD@iver
112 / 112 / 13
Регистрация: 01.10.2008
Сообщений: 876
15.07.2010, 10:06 2
у меня была такая же проблема была сортировка массива пузырьком
так у меня на асм было в 2 раза дольше чем на С++
оказалось что дело в команде xch, оказывается она выполнялас дльше чем 3 команды mov
короче советую тебе поискать у себя тяжелые команды, и заменять их
1
vital792
2002 / 1274 / 60
Регистрация: 05.06.2010
Сообщений: 2,213
15.07.2010, 14:39  [ТС] 3
Спасибо, но тут вроде нет медленных команд, таких как например xchg - блокирующих шину. Скорее всего данные просто не помещаются в кеш и время тратится на подгрузку. Короче придется изменять сам алгоритм или хотя бы синусы/косинусы предварительно поместить в таблицу.
0
ht1515
шарпопочитатель
58 / 25 / 7
Регистрация: 31.01.2010
Сообщений: 1,005
22.03.2011, 14:25 4
Но я ожидал что использование sse3 увеличит производительность и функция будет работать хотя бы вдвое быстрее аналога с сопроцессором.
а почему именно в 2 раза быстрее должна быстрее работать? Кстати с++ он же тоже в sse всякие переводит без ведома пользователя, для увеличения скорости.

Добавлено через 5 минут
а если на ffwt библиотеке сделать тоже самое быстрее будет?
0
murderer
3321 / 1467 / 134
Регистрация: 06.10.2010
Сообщений: 3,230
22.03.2011, 14:47 5
Обязательно нужна двойная точность?
0
vital792
2002 / 1274 / 60
Регистрация: 05.06.2010
Сообщений: 2,213
22.03.2011, 14:48  [ТС] 6
Цитата Сообщение от ht1515 Посмотреть сообщение
а почему именно в 2 раза быстрее должна быстрее работать?
хотя бы в 2 раза)). Не почему) просто я надеялся.
Цитата Сообщение от ht1515 Посмотреть сообщение
а если на ffwt библиотеке сделать тоже самое быстрее будет?
там преобразование фурье проводится в 2 этапа. Первый - построение плана (fftw_plan()), второй вызов самой функции fftw_execute(). Если сравнить по скорости эти 2 функции и мою - будет примерно одинаково, но если нужно вычислять fft в цикле (а это почти всегда - длинный сигнал разбивают на фрагменты, сглаживают оконной функцией и преобразуют) то быстродействие значительно возрастает, т.к. план строить достаточно один раз, а fft_execute() выполнять в цикле.
0
ht1515
шарпопочитатель
58 / 25 / 7
Регистрация: 31.01.2010
Сообщений: 1,005
22.03.2011, 14:50 7
murderer, какая разница ? Хоть тройная, он же алгоритмы сравнивает.
0
vital792
2002 / 1274 / 60
Регистрация: 05.06.2010
Сообщений: 2,213
22.03.2011, 14:55  [ТС] 8
Цитата Сообщение от murderer Посмотреть сообщение
Обязательно нужна двойная точность?
обязательно. Обычного float не хватает точности - при длинных сигналах накапливается ошибка округления. double дает вполне приемлемую погрешность
0
ht1515
шарпопочитатель
58 / 25 / 7
Регистрация: 31.01.2010
Сообщений: 1,005
22.03.2011, 14:56 9
vital792, то что у тебя на 15 % быстрее заработало, значит как компилятор ты круче того на котором писал вариант для с++.
0
vital792
2002 / 1274 / 60
Регистрация: 05.06.2010
Сообщений: 2,213
22.03.2011, 15:04  [ТС] 10
Цитата Сообщение от ht1515 Посмотреть сообщение
vital792, то что у тебя на 15 % быстрее заработало, значит как компилятор ты круче того на котором писал вариант для с++.
ht1515, этого недостаточно. Это даже и близко не сравнимо по скорости с fftw. А я надеялся именно ее заменить когда писал свою функцию. Оптимизировать надо сначала на уровне алгоритма, а потом в последнюю очередь на низком. Есть всякие методы (например вайнограда-фурье) которые очень сложны в реализации, но дают некоторый выигрыш в быстродействии, видимо в fftw учтено абсолютно все что можно - видно люди немало поработали. MTI forever!
0
ht1515
шарпопочитатель
58 / 25 / 7
Регистрация: 31.01.2010
Сообщений: 1,005
22.03.2011, 15:10 11
vital792, mit наверно
0
vital792
22.03.2011, 15:13  [ТС]
  #12

Не по теме:


Цитата Сообщение от ht1515 Посмотреть сообщение
mit наверно
да, извиняюсь. Пальцы кривые))

0
ISergey
Maniac
Эксперт С++
1413 / 923 / 149
Регистрация: 02.01.2009
Сообщений: 2,754
Записей в блоге: 1
22.03.2011, 22:29 13
vital792, Попробуй вычисления делать на видеокарте..
Здесь результаты вродь не плохие..
http://www.cv.nrao.edu/~pdemores/gpu/
1
ht1515
шарпопочитатель
58 / 25 / 7
Регистрация: 31.01.2010
Сообщений: 1,005
22.03.2011, 22:31 14
на гпу вроде какие то специализированные задачи решаю для графики, а не для математики... хотя хз.
0
murderer
3321 / 1467 / 134
Регистрация: 06.10.2010
Сообщений: 3,230
23.03.2011, 05:49 15
C
1
2
3
4
5
 for (bitm = 2, j = 0; bitm < 2*fftFrameSize; bitm <<= 1)
                {
                        if (i & bitm) j++;
                        j <<= 1;
                }
http://www.cyberforum.ru/assembler/t...ead233439.html
1
vital792
2002 / 1274 / 60
Регистрация: 05.06.2010
Сообщений: 2,213
23.03.2011, 08:43  [ТС] 16
Цитата Сообщение от ISergey Посмотреть сообщение
Попробуй вычисления делать на видеокарте..
да, я уже смотрел в сторону cuda, но к сожалению у меня не nvidia.

murderer, да, действительно оптимизировать немного еще можно, но основное время тратится не на этот цикл, а на следующий... В общем всем спасибо за советы, тема выплыла конечно поздновато, НИР уже сдана, но на будущее может пригодится.
0
ht1515
шарпопочитатель
58 / 25 / 7
Регистрация: 31.01.2010
Сообщений: 1,005
06.06.2011, 11:22 17
а пример вызова этой функции на с++ скажите как выглядит.

Хочу поэкспериментировать.
0
vital792
2002 / 1274 / 60
Регистрация: 05.06.2010
Сообщений: 2,213
06.06.2011, 11:49  [ТС] 18
вот так я сравнивал скорости fftw, fft на си и fft на асме:

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
#include "fftw3.h"
#include "math.h"
#include "windows.h"
 
#pragma comment(lib,"libfftw3-3.lib")
 
typedef double fftw_complex[2];
double pi = 3.1415926535897932384626433832795;
 
int size_data = 1 << 21;    //2^20
 
void Fft2(double *in_Buffer, double *out_Buffer, long fftFrameSize, long sign);
 
extern "C" void FFT_SSE2( double *in_data, double *out_data, int fftFrameSize);
 
int main(int argc, char* argv[])
{/*
    __declspec(align(16)) double *complex_s = (double*)_aligned_malloc(sizeof(double)*16, 16);
    double *_ptr = complex_s;
    for(int i=0; i<8; i++)
    {
        *(_ptr++) = (double)(i+1)/10;
        *(_ptr++) = (double)i+1.0;
    }
 
    __asm
    {
        mov     eax, complex_s;
        movapd  xmm0, xmmword ptr[eax+16]
        movddup xmm1, qword ptr [eax+8]
        movddup xmm2, qword ptr [eax]
        mulpd   xmm2, xmm0
        shufpd  xmm0, xmm0, 0x1
        mulpd   xmm1, xmm0
        addsubpd xmm2, xmm1
    }
*/
 
    fftw_plan p;
    double *signal = new double[size_data];// = {7,8,11,15,17,14,9,5,2,-1,-4};
    fftw_complex *out_signal = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * (size_data));// / 2 + 1));
 
    for(int i=0;i<size_data;i++)
        signal[i] = sin(2*pi*12800*i/250000);
 
 
double *smb_signal = (double*)_aligned_malloc(sizeof(double)*(size_data<<1), 16);//new double[size_data<<1];
double *temp = smb_signal;
 
for(int i=0; i<size_data; i++)
{
    *(temp++) = signal[i];
    *(temp++) = 0;
}
 
DWORD _time = GetTickCount();
    p = fftw_plan_dft_1d(size_data, (fftw_complex*)smb_signal, out_signal, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute(p);
_time = GetTickCount() - _time;
 
DWORD time_c = GetTickCount();
Fft2(smb_signal, (double*)out_signal, size_data, -1);
time_c = GetTickCount() - time_c;
 
DWORD time_asm = GetTickCount();
FFT_SSE2(smb_signal, (double*)out_signal, size_data);
time_asm = GetTickCount() - time_asm;
 
_aligned_free(smb_signal);
//delete [] smb_signal;
 
    fftw_destroy_plan(p);
    fftw_free(out_signal);
    delete [] signal;
    return 0;
}
 
void Fft2(double *in_Buffer, double *out_Buffer, long fftFrameSize, long sign)
{
    double wr, wi, arg, *p1, *p2, temp;
    double tr, ti, ur, ui, *p1r, *p1i, *p2r, *p2i;
    long i, bitm, j, le, le2, k;
 
    memcpy(out_Buffer, in_Buffer, fftFrameSize*2*sizeof(double));
 
    for (i = 2; i < 2*fftFrameSize-2; i += 2)
    {
        for (bitm = 2, j = 0; bitm < 2*fftFrameSize; bitm <<= 1)
        {
            if (i & bitm) j++;
            j <<= 1;
        }
        if (i < j)
        {
            p1 = out_Buffer+i; p2 = out_Buffer+j;
            temp = *p1; *(p1++) = *p2;
            *(p2++) = temp; temp = *p1;
            *p1 = *p2; *p2 = temp;
        }
    }
    long lim = (long)(log((double)fftFrameSize)/log(2.)+.5);
    for (k = 0, le = 2; k < lim; k++)
    {
        le <<= 1;
        le2 = le>>1;
        ur = 1.0;
        ui = 0.0;
        arg = pi / (le2>>1);
        wr = cos(arg);
        wi = sign*sin(arg);
        for (j = 0; j < le2; j += 2)
        {
            p1r = out_Buffer+j; p1i = p1r+1;
            p2r = p1r+le2; p2i = p2r+1;
            for (i = j; i < 2*fftFrameSize; i += le)
            {
                tr = *p2r * ur - *p2i * ui;
                ti = *p2r * ui + *p2i * ur;
                *p2r = *p1r - tr; *p2i = *p1i - ti;
                *p1r += tr; *p1i += ti;
                p1r += le; p1i += le;
                p2r += le; p2i += le;
            }
            tr = ur*wr - ui*wi;
            ui = ur*wi + ui*wr;
            ur = tr;
        }
    }
}
0
ht1515
шарпопочитатель
58 / 25 / 7
Регистрация: 31.01.2010
Сообщений: 1,005
06.06.2011, 12:12 19
и ещё не понятно что на выходе функции будет.
Спектр что ли?

Добавлено через 21 минуту
получается
C++
1
double *smb_signal = (double*)_aligned_malloc(sizeof(double)*(size_data<<1), 16);//new double[size_data<<1];
хмм чето как-то сложно для понимания.
если у меня есть массив исходных значений double mas[i] например. Я могу так сделать же?

double*smb_signal = mas[];



Что за тип такой? Заменить его С++шным можно?
C++
1
fftw_complex *out_signal = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * (size_data));// / 2 + 1));
Можно как предыдущую её объявить?

вызов
C++
1
Fft2(smb_signal, (double*)out_signal, size_data, -1);
а минус 1 что означает?
0
vital792
2002 / 1274 / 60
Регистрация: 05.06.2010
Сообщений: 2,213
06.06.2011, 12:40  [ТС] 20
Цитата Сообщение от ht1515 Посмотреть сообщение
и ещё не понятно что на выходе функции будет.
Спектр что ли?
На выходе будет фурье-образ сигнала. Из него можно получить спектр, например, если взять значения по модулю будет амплитудный спектр, если взять аргумент (arctg(im/re)) будет фазовый спектр.
Цитата Сообщение от ht1515 Посмотреть сообщение
хмм чето как-то сложно для понимания.
если у меня есть массив исходных значений double mas[i] например. Я могу так сделать же?
double*smb_signal = mas[];
Здесь выделяется динамическая память с выравниванием по 16 разрядной границе (команды sse при выравненной памяти работают быстрее).
Код
double*smb_signal = mas[]
так ты просто присвоишь указателю адрес своего массива, без всякого выравнивания памяти.

Цитата Сообщение от ht1515 Посмотреть сообщение
Что за тип такой? Заменить его С++шным можно?
Так ведь он определен в начале программы:
Цитата Сообщение от vital792 Посмотреть сообщение
typedef double fftw_complex[2];
Можно заменить просто структурой double[2]
Цитата Сообщение от ht1515 Посмотреть сообщение
а минус 1 что означает?
последний аргумент функции - тип преобразования, -1 прямое, 1 обратное
0
06.06.2011, 12:40
MoreAnswers
Эксперт
37091 / 29110 / 5898
Регистрация: 17.06.2006
Сообщений: 43,301
06.06.2011, 12:40

преобразование Фурье
Здравствуйте, Возник такой вопрос. На википедии дана формула Прямого...

Преобразование Фурье
Здравствуйте! Необходимо провернуть преобразование Фурье модельного сигнала...

Преобразование Фурье
Добрый день.Помогите. Нужна помощь нужно сделать преобразование Фурье. Нам дана...


Искать еще темы с ответами

Или воспользуйтесь поиском по форуму:
20
Ответ Создать тему
Опции темы

КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2018, vBulletin Solutions, Inc.
Рейтинг@Mail.ru