Wikipedia:Reference desk/Archives/Miscellaneous/2024 June 3
| 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
editSinc/Lanczos FFT bin interpolation on FFT-based log-frequency spectrum analyzers
editI 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 where 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)
- @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)
- 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)
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)