The transfer function can be calculated with any frequency accuracy, as opposed to a DFT or FFT.

There is a requirement that it must be a linear phase filter, ie with a symmetrical or antisymetric impulse response.

The algorithm is described in, among others, http://eeweb.poly.edu/iselesni/EL713/zoom/linphase.pdf

The code is in file 'DSP utils' that will contain multiple functions

Code: Select all

```
'DSP utils
'Functions for Digital Signal Processing
'by Ton Nillesen, 2019
'
DEF DSP_FIRgain(h(),f1,f2,df)
'by Ton Nillesen, October 2019
'Frequency-gain characteristic of linear phase FIR filter
'h(0)=filter coefficients, f1 to f2=freq range, df=freq step
'See http://eeweb.poly.edu/iselesni/EL713/zoom/linphase.pdf
'Error=10 if size of h()<2
'Error=11 if h() is not symmetric or anti-symmetric
'Error=12 if parameters are not accepted
'Function returns: List with gain data, max gain, min gain
N=SIZE OF h()
DIM A() ' frequency gain characteristic
IF N<2 : ERROR=10 : RETURN A() : ENDIF
IF (df<=0) OR (f2-f1)<=0 : ERROR=12 : RETURN A() : ENDIF
M=(N-1)/2
'--- check symmetry
FOR i=0 TO m
IF ABS(h(i))<>ABS(h(n-1-i)) : ERROR=11 : RETURN A()
ENDIF
NEXT i
'--- set type
GOSUB SetOdd
IF Odd
IF h(0)=h(n-1) THEN Type=1 ELSE Type=3
ELSE : IF h(0)=h(n-1) THEN Type=2 ELSE Type=4
ENDIF
'---- Presets
Amax=0 : Amin=0
IF (Type=1) or (Type=3) : imax=M-1
ELSE : imax=N/2-1
ENDIF
'--- process
FOR f=f1 TO f2 STEP df
A+0 : j=(SIZE OF A)-1 ' new entrie A(j)
IF type=1 THEN A(j)+=h(M) ' A=H(m) at f=0
IF Type<3 THEN GOSUB Cosine ELSE GOSUB Sine
NEXT f
RETURN A() Amax Amin ' from function
SetOdd: '--- local subroutine
IF INT(n/2)=n/2 THEN Odd=0 ELSE Odd=1
RETURN ' from subroutine
Cosine: '--- local subroutine
w=2PI*f ' angular frequency
FOR i=0 TO imax
A(j)+=2*h(i)*COS((M-i)*w)
Amax=MAX(Amax,A(j)) : Amin=MIN(Amin,A(j))
NEXT i
RETURN
Sine: '--- local subroutine
w=2PI*f ' angular frequency
FOR i=0 TO imax
A(j)+=2*h(i)*SIN((M-i)*w)
Amax=MAX(Amax,A(j)) : Amin=MIN(Amin,A(j))
NEXT i
RETURN
END DEF
```

Code: Select all

```
'DSP_FIRgain TEST, frequency-gain characteristic
'by Ton Nillesen, 2019
'See http://eeweb.poly.edu/iselesni/EL713/zoom/linphase.pdf
'
'==== Settings
f1=0 ' start frequency, relative to Fsample
f2=0.5 ' stop frequency, relative to Fsample
'==== Constants ====
SCRSIZE 600 400
sw sh = SCRSIZE|2
rb gb bb=0.8,1,0.8 ' background color
DRAW SIZE 1
radius=3
margin=50 ' left and right margin of plot
x0=margin ' start-x
y0=sh-2*margin ' baseline
ymax=0.9*y0 ' maximum y
'
'==== Examples of FIRtypes 1 to 4 ====
DIM h1()=(-0.047 ,0.101 ,-0.151 ,0.187 ,1.000 ,0.187 ,-0.151 ,0.101 ,-0.047)
DIM h2()=(0.042 ,-0.009 ,-0.050 ,0.053 ,0.034 ,-0.127 ,0.060 ,0.519 ,0.519 , 0.060 ,-0.127 ,0.034 ,0.053 ,-0.050 ,-0.009 ,0.042 )
DIM h3()=(-0.027,-0.055,-0.089,-0.133,-0.200,-0.322,-0.670,0.000,0.670,0.322,0.200,0.133,0.089,0.055,0.027)
DIM h4()=(-0.1,0.6,-0.6,0.1)
'
'==== Main ====
GOSUB Init
DRAW CLEAR rb,gb,bb
DRAW COLOR 1,0,0
'--- show all types
Type=1
DO IF Type<5
INSIGNTEXT 0 "FIR type "+STRTEXT type
DIM H()
IF type=1 THEN H=h1
IF type=2 THEN H=h2
IF type=3 THEN H=h3
IF type=4 THEN H=h4
DIM A()
A Amax Amin = DSP_FIRgain|3(H,f1,f2,(f2-f1)/(sw-2*margin))
n1=(SIZEOF a)-1
GOSUB Rectify
GOSUB Chars
GOSUB Plot
WaitForNext
Type+=1
REDO
DRAW CLEAR 1,1,1
INSIGNTEXT 0 "Done"
INSIGN POS 0 0,sh/2
INBUTT HIDE 0
END
'
'==== Functions and subroutines
{DSP utils}
{interface 1.part}
'
Plot:
DRAW CLEAR rb,gb,bb
dx=(sw-2*margin)/n1
dy=ymax/Amax
'--- axes
DRAW COLOR 0,0,0
DRAW SIZE 2
DRAW LINE x0,y0,x0+n1*dx,y0
DRAW LINE x0,y0,x0,0
'--- graduation
DRAW FONT "Menlo",16
th=STRHEIGHT "X"
'--- x-marks
x1=FLOOR(10*f1)/10
x2=CEIL(10*f2)/10
FOR i=x1 TO x2 STEP 0.1
txt$=(STRTEXT i "0.0")
tw=STRWIDTH txt$
x=(i-x1)*n1/(f2-f1)
DRAW TEXT txt$,x0+x-tw/2, y0+2*th
DRAW LINE x0+x,y0,x0+x,y0+th/2
DRAW SIZE 1
DRAW LINE x0+x,y0,x0+x,0
NEXT i
'--- y-marks
DRAW SIZE 2
atop=CEIL(10*Amax)/10
FOR i=0 TO atop STEP 0.2
txt$=(STRTEXT i "0.0")
tw=STRWIDTH txt$
y=y0-i*dy
DRAW TEXT txt$,x0-1.5*tw, y+th/2
DRAW LINE x0,y,x0-th/2,y
DRAW SIZE 1
DRAW LINE x0,y,x0+n1*dx,y
NEXT i
'--- curve
DRAW COLOR 1,0,0
DRAW SIZE 3
DRAW AT x0,y0-dy*A(0)
FOR i=0 TO n1
x=x0+i*dx : y=y0-dy*A(i)
DRAW TO x,y
NEXT i
RETURN
'
Rectify: ' make positive and set max value
FOR i=0 TO n1 : A(i)=ABS(A(i)) : NEXT i
IF n1<25 THEN 'print list
PRINT "List A size =";n1+1
FOR i=0 TO n1
PRINT A(i)
NEXT i
ENDIF
Amax=MAX(Amax,ABS(Amin)) ' rectified
Amin=0
RETURN
'
Chars: ' print characteristics
PRINT "FIR type="+(STRTEXT DSP_FIRgain.type "?, ");
PRINT " M="+(STRTEXT DSP_FIRgain.m "?, ");
PRINT " Size="+(STRTEXT DSP_FIRgain.n "? coëfficients")
PRINT "Error="+(STRTEXT ERROR)
PRINT "Frequency-gain list A =";n1+1;"points"
PRINT "Amax =";Amax : PRINT "Amin =";Amin
PRINT
RETURN
```

Attaching images gave problems, however MrK. solved that and posted the images.

See next post by kibernetik » Wed Oct 23, 2019 8:23 am