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

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

14.07.2010, 10:44. Просмотров 6068. Ответов 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
Я подобрал для вас темы с готовыми решениями и ответами на вопрос Преобразование Фурье на ассемблере (Assembler):

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

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

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

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

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

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

43
kosem
0 / 0 / 0
Регистрация: 07.12.2011
Сообщений: 11
08.12.2011, 14:32 #41
Спасибо Виталий, буду тестировать Matlab есть про степень двойки знаю.
0
kosem
0 / 0 / 0
Регистрация: 07.12.2011
Сообщений: 11
09.12.2011, 23:32 #42
Виталий, присланная Вами программа у меня тоже падает вот на этом месте. Возможно Вы знаете в чем дело если она работает на Вашей машине. Рис в низу.
0
Миниатюры
Преобразование Фурье на ассемблере  
vital792
1997 / 1269 / 60
Регистрация: 05.06.2010
Сообщений: 2,213
10.12.2011, 10:07  [ТС] #43
да, правда, извиняюсь, у меня тоже не работает. Функция 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
Программа на с++ с исправленной функцией во вложении
1
Вложения
Тип файла: zip tst_fft.zip (44.7 Кб, 35 просмотров)
kosem
0 / 0 / 0
Регистрация: 07.12.2011
Сообщений: 11
10.12.2011, 15:01 #44
Благодарю Вас Виталий, все работает нормально.Респект.
0
10.12.2011, 15:01
MoreAnswers
Эксперт
37091 / 29110 / 5898
Регистрация: 17.06.2006
Сообщений: 43,301
10.12.2011, 15:01
Привет! Вот еще темы с решениями:

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

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

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

Преобразование Фурье
Есть программа в которую нужно добавить Преобразование Фурье, и вывести...


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

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

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