Wikipedia:Reference desk/Archives/Miscellaneous/2024 June 3

= 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)
 * Isn't this more of a maths desk kind of thing? --Viennese Waltz 07:06, 4 June 2024 (UTC)
 * 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)
 * @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)
 * 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)

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)