gens
Joined: 18 Feb 2013
Posts: 161

FIR low pass filter
stereo, 16bit
44100 Hz sampling frequency
4000 Hz cutoff frequency
Blackman window
32 taps
fixed point
attuned to 0dB
forgot to add
amd64 and linux
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
frequency response, needs some work (but not much)
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*(nm)) ) / (pi*(nm))
hd[n]
}
#print "BlackmanHarris 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
filter calculator in bc language
needs bc ofc, run it with bc l filter_calc.bc
Last edited by gens on 28 Mar 2014, 02:58; edited 1 time in total
