format ELF64 executable
entry start

align 16
segment readable writeable

samples equ 1024

align 16
buff_in: rw samples*2
align 16
buff_out: rw samples*2

; 2^16 = 1.0
align 16
coefficients:

dw 5
dw 14
dw 35
dw 57
dw 39
dw -76
dw -328
dw -677
dw -945
dw -828
dw 16
dw 1799
dw 4428
dw 7441
dw 10109
dw 11680
dw 11680
dw 10109
dw 7441
dw 4428
dw 1799
dw 16
dw -828
dw -945
dw -677
dw -328
dw -76
dw 39
dw 57
dw 35
dw 14
dw 5

rq 10 ;just in case
tmp rq 100

sys_read equ 0
sys_write equ 1
sys_exit equ 60


align 16
segment readable executable

start:

xorps xmm7, xmm7
xorps xmm6, xmm6
xorps xmm5, xmm5
xorps xmm4, xmm4
xorps xmm3, xmm3
xorps xmm2, xmm2
xorps xmm1, xmm1
xorps xmm0, xmm0

mov [tmp], rdx
mov [tmp+8], rdx

filter:

mov rdx, samples*2*2
mov rsi, buff_in
mov rdi, 0
mov rax, sys_read
syscall
cmp rax, 0
jz end_fir

xor rax, rax

mov rcx, samples
loopy:
movaps xmm12, [coefficients]
movaps xmm13, [coefficients+16*1]
movaps xmm14, [coefficients+16*2]
movaps xmm15, [coefficients+16*3]
mov edx, [rax+buff_in] ;dx = left channel sample
mov esi, edx
shr esi, 16 ;si = right channel sample

; left channel delay line
pslldq xmm3, 2
movaps xmm10, xmm2
psrldq xmm10, 2*7
orpd xmm3, xmm10

pslldq xmm2, 2
movaps xmm10, xmm1
psrldq xmm10, 2*7
orpd xmm2, xmm10

pslldq xmm1, 2
movaps xmm10, xmm0
psrldq xmm10, 2*7
orpd xmm1, xmm10

pslldq xmm0, 2 ;bytes
movzx rdx, dx
movq xmm10, rdx
orpd xmm0, xmm10

; right channel delay line
pslldq xmm7, 2
movaps xmm10, xmm6
psrldq xmm10, 2*7
orpd xmm7, xmm10

pslldq xmm6, 2
movaps xmm10, xmm5
psrldq xmm10, 2*7
orpd xmm6, xmm10

pslldq xmm5, 2
movaps xmm10, xmm4
psrldq xmm10, 2*7
orpd xmm5, xmm10

pslldq xmm4, 2
movq xmm10, rsi
orpd xmm4, xmm10




movaps xmm8, xmm12
movaps xmm9, xmm13
movaps xmm10, xmm14
movaps xmm11, xmm15


pmaddwd xmm12, xmm0
pmaddwd xmm13, xmm1
pmaddwd xmm14, xmm2
pmaddwd xmm15, xmm3

paddd xmm12, xmm14
paddd xmm13, xmm15

paddd xmm12, xmm13

movhlps xmm13, xmm12
paddd xmm12, xmm13

movss xmm13, xmm12
psrldq xmm12, 32
paddd xmm12, xmm13



pmaddwd xmm8, xmm4
pmaddwd xmm9, xmm5
pmaddwd xmm10, xmm6
pmaddwd xmm11, xmm7

paddd xmm8, xmm10
paddd xmm9, xmm11

paddd xmm8, xmm9

movhlps xmm9, xmm8
paddd xmm8, xmm9

movss xmm9, xmm8
psrldq xmm8, 32
paddd xmm8, xmm9



movd r10d, xmm12
;sal r10d, 3 ;overflow possibility unless compensated
shr r10d, 16
mov [rax+buff_out], r10w

movq r10, xmm8
;sal r10d, 2
shr r10d, 16
mov [rax+buff_out+2], r10w

add rax, 4
dec rcx
jnz loopy

mov rdx, samples*2*2
mov rsi, buff_out
mov rdi, 1
mov rax, sys_write
syscall

jmp filter

end_fir:
mov rax, sys_exit
syscall    

define int(x) {
    auto old_scale   /* variables global by default */
    old_scale=scale  /* scale is global */
    scale=0; ret=x/1
    scale=old_scale
    return ret
}
define round(x) {
    if (x<0) x-=.5 else x+=.5
    return int(x)
}

scale=50

#pi=3.14159265358979323846264338327950288
pi=4*a(1)
print "pi is "
pi
print "\n"

fs = 44100
fc = 4000
omega = (2*pi*fc)/fs
omega/pi

nf=31
nt=nf+1
m=nf/2

print "ideal filter coefficients\n"
for (n=0; n<nt; n++) {
    if (n == m) {
    hd[n]=(omega/pi)
    hd[n]
    continue }
    
    hd[n]=( s(omega*(n-m)) ) / (pi*(n-m))
    hd[n]
}


#print "Blackman-Harris window coefficients\n"
#for (n=0; n<nt; n++) {
#    w[n] = 0.35875 - 0.48829*c( (2*pi*n)/nf ) + 0.14128*c( (4*pi*n)/nf ) - 0.01168*c( (6*pi*n)/nf )
#    w[n]
#}

print "Blackman window coefficients\n"
for (n=0; n<nt; n++) {
    w[n] = 7938/18608 -  (9240/18608)*c( (2*n*pi)/nf ) + (1430/18608)*c( (4*n*pi)/nf )
    w[n]
}

#print "Blackman window coefficients\n"
#for (n=0; n<nt; n++) {
#    w[n] = 0.42 -  0.5*c( (2*n*pi)/nf ) + 0.08*c( (4*n*pi)/nf )
#    w[n]
#}

for (n=0; n<nt; n++) {
    h[n] = w[n]*hd[n]
    h[n]
}

print "coefficients sum\n"
for (n=0; n<nt; n++) {
    coeff_sum+=h[n]
}
coeff_sum
print "one over coefficients sum\n"
1/coeff_sum

print "final filter coefficients\n"
for (n=0; n<nt; n++) {
    h[n] = h[n] / coeff_sum
    #h[n]
}

print "normalised to 2^16\n"
for (n=0; n<nt; n++) {
    hn[n] = round( h[n]*(2^16) )
    print "dw "
    hn[n]
    total = total + hn[n]
}

print "new sum\n"
total
print "aim\n"
2^16