Miscellaneous desk
< June 2	<< May \| June \| Jul >>	Current desk >

Welcome to the Wikipedia Miscellaneous Reference Desk Archives
The page you are currently viewing is a transcluded archive page. While you can leave answers for any questions shown below, please ask new questions on one of the current reference desk pages.

June 3

Sinc/Lanczos FFT bin interpolation on FFT-based log-frequency spectrum analyzers

I wonder if sinc/Lanczos interpolation is the best FFT bin interpolation method for "bandpower" spectrum (especially on lower frequencies part when using logarithmic frequency scale) because it approximates $\int _{l}^{h}|{\widehat {f}}(x)|\,d\omega$ where ${\widehat {f}}(x)$ is a discrete-time Fourier transform, in the best way when summation mode is set to "Sum" on my relevant CodePen project or is it? 2001:448A:3070:E3DA:7021:FEBA:971D:F9F6 (talk) 22:27, 3 June 2024 (UTC)[reply]

Isn't this more of a maths desk kind of thing? --Viennese Waltz 07:06, 4 June 2024 (UTC)[reply]

Not really. It is a question that audio engineers might be able to answer if they are knowledgeable about digital signal processing. In this context, there is no mathematical notion of how "good" a technique is. I know what FFT is, I know what interpolation is, but not what "FFT bin interpolation" is. You will not find the term "FFT bin" in a maths handbook. --Lambiam 14:46, 4 June 2024 (UTC)[reply]

@Lambiam What I meant by "FFT bin interpolation" is the same interpolation is done on FFT bins, which is necessary for logarithmic frequency spectrum analyzers since FFT have limited resolution on lower frequencies (since it has linear frequency resolution as opposed to frequency bands, which in this case has logarithmic frequency scale) and sinc interpolation closely approximates the zero-padding I believe when the interpolation is done before conversion to magnitude FFT. 2001:448A:3070:E3DA:E523:AA53:7234:600A (talk) 23:47, 4 June 2024 (UTC)[reply]

Yes. Wikipedia uses extended meanings of the useful term BIN for a partition or discrete interval in a range of values such as a Histogram bin, Data binning, a data pre-processing technique or Bin (computational geometry) a space partitioning data structure to enable fast region queries and nearest neighbor search. Philvoids (talk) 22:31, 5 June 2024 (UTC)[reply]

It's unlikely that a volunteer here will shepherd the OP in their coding project on another web site but I can give references and a worked example that will be useful. The IEC standard 61672 for sound level meters gives much information on professional-standard audio spectral analysis.

A spectrum analysis in 1/3-octave steps implies using this bank of filters:

FILTER FREQUENCIES in Hz
       Band limits
Center Lower Upper

100   89.1   112
125    112   141
160    141   178
200    178   224
250    224   282
315    282   355
400    355   447
500    447   562
630    562   708
800    708   891
1000   891  1122
1250  1122  1413
1600  1413  1778
2000  1778  2239
2500  2239  2818
3150  2818  3548
4000  3548  4467
5000  4467  5623
6300  5623  7079
8000  7079  8913
10000 8913 11220

I distinguish two alternative approaches. 1) Bank of disparate filters, or 2) Single FFT with tailored bin allocations. In either case the effort in the project depends on the chosen goals for resolution in power and frequency, and whether a near real-time spectrum display is required. (Performing the latter with high precision demands dedicated hardware such as a DSP or FPGA.)

1) Bank of disparate filters The table defines 21 bandpass filters each of different widths, to be separately designed. The Butterworth bandpass design is optimised for flat response between its lower and upper half-power (-3dB) points. This means that other filter characteristics such as the out-of-band rolloff rates are neglected or "poorly" shaped. The best we can do is let adjacent filters overlap at their -3dB frequencies. Here is a worked example in Fortran to design a single Butterworth bandpass filter.

=====================================================

Example requirements

Change as required for each filter.
Band: 20 to 30 kHz
Sections: 5
Sampling interval: 10 usec

	DIMENSION A(5),B(5),C(5),D(5),E(5),GRAF(2.20)
	CALL BPDES(20000.,30000.,1.E-5,5,A,B,C,D,E,GRAF)
	DO 1 N=1,5
1	WRITE(5,2) N,A(N),B(N),C(N),D(N),E(N)
2	FORMAT(5(13,5E14.6))
	DO 3 N=1,20
	DB=10.*ALOG10(GRAF(2,N))
3	WRITE(5,4) GRAF(1,N),GRAF(2,N),DB
4	FORMAT(1X,3E15.6)
	STOP
	END

C - BPDES
C - BANDPASS BUTTERWORTH DIGITAL FILTER DESIGN SUBROUTINE
C - INPUTS ARE PASSBAND (3-DB) FREQUENCIES F1 AND F2 IN HZ.
C -            SAMPLING INTERVAL T IN SECONDS, AND
C -            NUMBER NS OF FILTER SECTIONS.
C - OUTPUTS ARE NS SETS OF FILTER COEFFICIENTS, I.E.
C -                 A(K) THRU E(K) FOR K=1 THRU NS, AND
C-              20 PAIRS OF FREQUENCY AND POWER GAIN, I.E.
C -                 GRAF(1,K) aND GRAF(2,K) FOR K=1 THRU 20.
C - NOTE THAT A(K) THRU E(K) AS WELL AS GRAF(2,20) MUST BE
C -      DIMENSIONED IN THE CALLING PROGRAM.
C -
C - THE DIGITAL FILTER HAS NS SECTIONS IN CASCADE. THE KTH
C -     SECTION HAS THE TRANSFER FUNCTION
C -
C -                  A(K)+(Z**4-2*Z**2+1)
C -     H(Z)=--------------------------------------
C -           Z**4+B(K)*Z**3+C(K)*Z**2+D(K)*Z+E(K)
C -
C - THUS, IF F(M) and G(M) ARE THE INPUT AND OUTPUT OF THE
C -     KTH SECTION AT TIME M*T, THEN
C - 
C -   G(M)=A(K)*(F(M)-2*F(M-2)+F(M-4))-B(K)*G(M-1)
C -          -C(K)*G(M-2)-D(K)*G(M-3)-E(K)*G(M-4)
C -
	SUBROUTINE BPDES(F1,F2,T,NS,A,B,C,D,E,GRAF)
	DIMENSION A(1),B(1),C(1),D(1),E(1),GRAF(2,20)
	PI=3.1415926536
	W1=SIN(F1*PI*T)/COS(F1*PI*T)
	W2=SIN(F2*PI*T)/COS(F2*PI*T)
	WC=W2-W1
	Q=WC*WC+2.*W1*W2
	S=W1*W1*W2*W2
	DO 150 K=1,NS
	CS=COS(FLOAT(2*(K+NS)-1)*PI/FLOAT(4*NS))
	P=-2.*WC*CS
	R=P*W1*W2
	X=1.+P+Q+R+S
	A(K)=WC*WC/X
	B(K)=(-4.-2.*P+2.*R+4.*S)/X
	C(K)=(6.-2.*Q+6.*S)/X
	D(K)=(-4.+2.*P-2.*R+4.*S)/X
150	E(K)=(1.-P+Q-R+S)/X
	DO 160 J=1,2
	DO 160 I=1,10
	K=I*(2-J)+(21-I)*(J-1)
	GRAF(2,K)=.01+.98*FLOAT(I-1)/9
	X=(1./GRAF(2,K)-1.)**(1./FLOAT(4+NS))
	SQ=SQRT(WC*WC*X*X+4.*W1*W2)
160	GRAF(1,K)=ABS(ATAN(.5*(WC*X+FLOAT(2*J-3)*SQ)))/(PI*T)
	RETURN
	END

=====================================================

The first WRITE statement above lists the coefficients of the 5-section filter:

K	A(K)	    	B(K)	C(K)		D(K)	E(K)

1	0.87451E-01	-0.*	0.14818E+01	-0.*	0.83158E+00
2	0.75377E-01 	-0.*	0.12772E+01	-0.*	0.57872E+00
3	0.67455E-01	-0.*	0.11430E+01	-0.*	0.41280E+00
4	0.62671E-01	-0.*	0.10619E+01	-0.*	0.31258E+00
5	0.60417E-01	-0.*	0.10237E+01	-0.*	0.26538E+00

Values of B(K) and D(K) are theoretically zero in this case but show tiny rounding errors.

The second WRITE statement lists 20 points on the power gain curve of the resulting filter.

This confirms -3dB gains at F1 and F2 and you can assess the overlap between adjacent filters.

FREQ (HZ)	       POWER GAIN	POWER GAIN (DB)
 	
10 points on lower	::::::::	::::::::
skirt including F1 

10 points on upper
skirt including F2	::::::::	::::::::
=====================================================

The reference for the above program is "Digital Signal Analysis" by Samuel D. Stearns, 1975 Hayden. An updated version that includes an IBM floppy disc is "Digital Signal Analysis" by Samuel D. Stearns and Don R. Rush, 1990 Hayden.

2) Single FFT with tailored bin allocations

A possible FFT specification:

	Sample rate: 32,768 samples/second
	Inputs: 3,276 real, imaginary components zero
	Outputs: 3,276 power squared (I² + Q²)
This implements 3,275 bandpass filters 10 times a second.

To obtain the overall level in any bandwidth, sum the squared levels of each FFT bin in the band, divide by the number of bins, and take the square root. Where the FFT bin -3dB points do not match the desired 1/3-octave steps, share the bin powers across borders as in this example:

FREQUENCIES in Hz
1/3-OCTAVE FILTER (example)   
       Band limits           ADD THESE FFT BINS   Weight
Center Lower Upper           Center Lower  Upper
125    112                   110    105    115     30%
                             120    115    125    100%
                             130    125    135    100%
             141             140    135    145     40%

More analysis filter resolution at low frequencies will need FFT process of longer sound streams than 0.1 second while analysis to higher audio frequencies will need a higher sample rate. Pursuing both these aims will increase demand on the FFT computation, give a slower result or need more expensive hardware. I do not think that alternative filter designs from the analog world offer a shortcut. Philvoids (talk) 23:25, 5 June 2024 (UTC)[reply]