# Universal filtering
# Made by AudiTon (auditon@xs4all.nl) october 2007
# 1st order high pass and 1st order low pass sections in cascade
# time constants can be filled-in
# result: 1. selected Sound filtered
#         2. filter available as spectrum object (complex)
#         3. function can be drawn if asked for
#         4. filter can be made by analog hardware

form Universal filter, made by 1st order lp and hp sections
    comment (Disable filter section by filling-in zero time constant)
    comment 
    real Time_constant_low_pass_1_(sec) 5e-4
    real Time_constant_low_pass_2_(sec) 5e-4
    real Time_constant_low_pass_3_(sec) 5e-4
    real Time_constant_low_pass_4_(sec) 5e-4
    real Time_constant_high_pass_1_(sec) 0
    real Time_constant_high_pass_2_(sec) 0
    real Time_constant_high_pass_3_(sec) 0
    real Time_constant_high_pass_4_(sec) 0
    boolean Plot_filter_gain_function yes
endform

t1 = time_constant_low_pass_1
t2 = time_constant_low_pass_2
t3 = time_constant_low_pass_3
t4 = time_constant_low_pass_4
t5 = time_constant_high_pass_1
t6 = time_constant_high_pass_2
t7 = time_constant_high_pass_3
t8 = time_constant_high_pass_4

# make base spectrum
obj$ = selected$("Sound")
To Spectrum... yes
Copy... former
# make modulus 1
Formula... 1/sqrt(2)

# low pass filters
for i to 4
    if t'i' != 0
        tc= t'i'
        call lowpass
    endif
endfor

# high pass filters

for i from 5 to 8
    if t'i' != 0
        tc= t'i'
        call highpass
    endif
endfor

Copy... filter

# filter sound
re$="Spectrum_'obj$'[1,col]"
im$="Spectrum_'obj$'[2,col]"
select Spectrum former
call complexmult
Rename... filtered
To Sound

if plot_filter_gain_function=1
    select Spectrum filter
    # make unity gain ref. 0 dB in picture
    # When modulus = 1 then level = 94 dB SPL
    bin=Get bin width
    Formula... self*10^(-94/20)/sqrt(2*'bin')
    Draw (log freq)... 10 20000 -60 6 yes
    One mark left... 0 yes yes yes  
    # adjust level back to normal for possible more filtering
    Formula... self*10^(94/20)*sqrt(2*'bin')
endif



procedure lowpass
    re$="1/(1+(2*pi*x)^2*'tc'^2)"
    im$="-2*pi*x*'tc'/(1+(2*pi*x)^2*'tc'^2)"
    call complexmult
endproc

procedure highpass
    re$="(2*pi*x)^2*'tc'^2/(1+(2*pi*x)^2*'tc'^2)"
    im$="2*pi*x*'tc'/(1+(2*pi*x)^2*'tc'^2)"
    call complexmult
endproc

procedure complexmult
    Copy... filter
    Formula... if row=1 then Spectrum_former[1,col]*'re$'-Spectrum_former[2,col]*'im$'
    ... else Spectrum_former[1,col]*'im$'+Spectrum_former[2,col]*'re$' fi
    select Spectrum former
    Remove
    select Spectrum filter
    Rename... former
endproc