## FIR gain

For finished programs
Dutchman
Posts: 151
Joined: Tue Aug 06, 2019 4:47 pm
Location: Netherlands

### FIR gain

I have created a DSP function that calculates the frequency transfer of a FIR filter.
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
``````
The test program is:

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
``````
The testprogram includes "interface 1.txt" which is attached
Interface 1.txt
(920 Bytes) Downloaded 78 times
.
Attaching images gave problems, however MrK. solved that and posted the images.
See next post by kibernetik » Wed Oct 23, 2019 8:23 am
Last edited by Dutchman on Thu Oct 24, 2019 3:08 pm, edited 3 times in total.
It is still a long way to go

Dutchman
Posts: 151
Joined: Tue Aug 06, 2019 4:47 pm
Location: Netherlands

### Re: FIR gain

See next post for images
Last edited by Dutchman on Thu Oct 24, 2019 3:07 pm, edited 1 time in total.
It is still a long way to go

kibernetik
Site Admin
Posts: 147
Joined: Tue Aug 06, 2019 3:03 pm

### Re: FIR gain

Very interesting program!

I did not understand what was wrong with screenshots posting.
Here are all 4 screenshots exactly as they were exported by screenshot operation (screenshot DPI setting=72): 1.png (63.63 KiB) Viewed 1943 times 2.png (66.96 KiB) Viewed 1943 times 3.png (73.5 KiB) Viewed 1943 times 4.png (66.58 KiB) Viewed 1943 times

Dutchman
Posts: 151
Joined: Tue Aug 06, 2019 4:47 pm
Location: Netherlands

### Re: FIR gain

While making the post, including two code listings, I could not add more after the first 3 attachments: 1 text and 2 pictures.
So all in all it probably went to a limit and was more than in your post It is still a long way to go

kibernetik
Site Admin
Posts: 147
Joined: Tue Aug 06, 2019 3:03 pm

### Re: FIR gain

Dutchman wrote:
Wed Oct 23, 2019 8:18 am
While making the post, including two code listings, I could not add more after the first 3 attachments: 1 text and 2 pictures.
So all in all it probably went to a limit and was more than in your post Ok, I see.
It appeared to be an attachment limit and it is increased to 15 items.

Dutchman
Posts: 151
Joined: Tue Aug 06, 2019 4:47 pm
Location: Netherlands

### Re: FIR gain

Function DSP_FIRgain has been simplified and its file name changed to "DSP utils"
Function now delivers the list with gain data and two additional values with maximum and minimum gain.
The test program has been adapted.
It is still a long way to go