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

@vital792 · Регистрация: 05.06.2010

Студворк — интернет-сервис помощи студентам

Ну так и здравствуйте! Имеется код на с (функция выполняет быстрое преобразование фурье):

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

@kosem · 08.12.2011, 14:32

Спасибо Виталий, буду тестировать Matlab есть про степень двойки знаю.

@kosem · 09.12.2011, 23:32

Виталий, присланная Вами программа у меня тоже падает вот на этом месте. Возможно Вы знаете в чем дело если она работает на Вашей машине. Рис в низу.

@vital792 · 10.12.2011, 10:07 **[ТС]**

да, правда, извиняюсь, у меня тоже не работает. Функция FFT_SSE2 немного не корректна: на выходе портятся регистры edi и esi. В debug версии на с++ это не играет роли, но в release получается что адреса массивов все время работы лежат в этих регистрах. Код функции с исправлением:

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
166
167
168
169
170
.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
sub     esp, 28
push    ebx
push    esi
push    edi
 
; 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
 
pop edi
pop esi
pop ebx
mov     esp, ebp
pop     ebp
retn
FFT_SSE2 endp
 
end

Программа на с++ с исправленной функцией во вложении

@kosem · 10.12.2011, 15:01

Благодарю Вас Виталий, все работает нормально.Респект.

Новые блоги и статьи Все статьи Все блоги /
Nekobox - outbounds[0].transport: unknown transport type: raw damix 01.10.2026 Фикс ошибки Правым кликом по серверу -> отладочная информация -> edit Заменить "net": "raw", на "net": "tcp", Нажать кнопку reload.	Программный домашний кинотеатр russiannick 27.09.2026 Сподобился на программный домашний кинотеатр. В качестве ЯВУ по традиции выбрал js. В помощники взял Яндекс-Алису. Было создано три зала на разные интересы. исторические и ретро сериал Хичкок. . .	Беседа с ИИ о программистах, недопускающих к созданию и правке кода генеративные ИИ и причины этого zorxor 21.09.2026 Раньше я радовался или получал некоторые эмоции, пусть небольшие, но всё же, от самого процесса написания кода, рекомпиляции и запуска, видя постепенное развитие программы и прочее. А теперь лень. . .	Мобильное приложение ColorStep pavlinmavlin 17.09.2026 Реализовал приложение Красный, Зеленый, Синий в Unity3d + c#. Название изменил на ColorStep. Приложение прошло модерацию и теперь доступно для скачивания. Делал его сам, шаг за шагом — и вот,. . .
Запрет дублирования строк в табличной части Maks 13.09.2026 Реализация из решения ниже выполнена на нетиповом справочнике "Нормы ТО" с табличной часть "Виды ТО", разработанного в КА2, со следующими реквизитами: - ВидТО (СправочникСсылка. ВидыТО); - ВидГСМ. . .	Скрипты Tampermonkey для CyberForum, ChatGPT, Claude и пр. Jin X 06.09.2026 Скрипты Tampermonkey для CyberForum, ChatGPT, Claude и пр. Работая с форумом и нейросетями в браузере часто хочется что-то подкорректировать или добавить какого-то функционала. Ниже прикреплён. . .	Программа опроса у.з. расходомера SLS-720F Argus19 02.09.2026 Программа опроса у. з. расходомера SLS-720F Программа опрашивает один раз в минуту три ультразвуковых расходомера SLS-720F через интерфейс RS-485 по протоколу Modbus RTU. Опрашиваются регистры. . .	Hyper-V: Компьютер должен поддерживать доверенный платформенный модуль 2.0. Maks 31.08.2026 При установке Windows 11 на виртуальную машину Hyper-V 2-го поколения вылезла такая ошибка: Решение: в параметрах виртуальной машины, в разделе "Безопасность" (Security) активировать флаг. . .

@kosem 0 / 0 / 0 Регистрация: 07.12.2011 Сообщений: 11
	08.12.2011, 14:32
	Спасибо Виталий, буду тестировать Matlab есть про степень двойки знаю. 0

@kosem 0 / 0 / 0 Регистрация: 07.12.2011 Сообщений: 11
	09.12.2011, 23:32
	Виталий, присланная Вами программа у меня тоже падает вот на этом месте. Возможно Вы знаете в чем дело если она работает на Вашей машине. Рис в низу. Миниатюры 0

@kosem 0 / 0 / 0 Регистрация: 07.12.2011 Сообщений: 11
	10.12.2011, 15:01
	Благодарю Вас Виталий, все работает нормально.Респект. 0

Опции темы