FIR gain

For finished programs
Post Reply
User avatar
Dutchman
Posts: 148
Joined: Tue Aug 06, 2019 4:47 pm
Location: Netherlands

FIR gain

Post by Dutchman » Tue Oct 22, 2019 7:33 pm

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 7 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

User avatar
Dutchman
Posts: 148
Joined: Tue Aug 06, 2019 4:47 pm
Location: Netherlands

Re: FIR gain

Post by Dutchman » Tue Oct 22, 2019 7:44 pm

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

User avatar
kibernetik
Site Admin
Posts: 131
Joined: Tue Aug 06, 2019 3:03 pm

Re: FIR gain

Post by kibernetik » Wed Oct 23, 2019 6:23 am

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
1.png (63.63 KiB) Viewed 257 times
2.png
2.png (66.96 KiB) Viewed 257 times
3.png
3.png (73.5 KiB) Viewed 257 times
4.png
4.png (66.58 KiB) Viewed 257 times

User avatar
Dutchman
Posts: 148
Joined: Tue Aug 06, 2019 4:47 pm
Location: Netherlands

Re: FIR gain

Post by Dutchman » 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 ;)
It is still a long way to go

User avatar
kibernetik
Site Admin
Posts: 131
Joined: Tue Aug 06, 2019 3:03 pm

Re: FIR gain

Post by kibernetik » Wed Oct 23, 2019 8:39 am

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.

User avatar
Dutchman
Posts: 148
Joined: Tue Aug 06, 2019 4:47 pm
Location: Netherlands

Re: FIR gain

Post by Dutchman » Thu Oct 24, 2019 3:11 pm

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

Post Reply