Blame src/matlab/wv_peakfreq.m at main - tutorials/wave-acoustics

2% WV_PEAKFREQ The frequency of the strongest partial of a between flo and fhi.

4% This is how every frequency in the tutorial is measured, from a recorded

5% signal and nothing else: taper the record with wv_window, find the largest

6% peak of its zero-padded DFT inside the band, and then polish the peak with

7% Newton's method on |X(f)|^2, using the exact derivatives of the windowed

8% discrete-time Fourier transform. For a clean partial well separated from

9% its neighbours the result is good to about 1e-10 relative; the method used

10% on it never needs to be.

12% a is a real record sampled every dt seconds; the band is in hertz.

13a = a(:);

14n = numel(a);

15w = wv_window(n);

16b = w .* (a - mean(a));

17nfft = 2 ^ nextpow2(4 * n);

18B = abs(fft(b, nfft));

19df = 1 / (nfft * dt);

20k0 = max(2, floor(flo / df) + 1);

21k1 = min(floor(nfft / 2), ceil(fhi / df) + 1);

22[~, i] = max(B(k0:k1));

23f = (k0 + i - 2) * df;

24tt = (0:n-1)' * dt;

25for it = 1:30

26 e = exp(-2i * pi * f * tt);

27 X = sum(b .* e);

28 X1 = sum(b .* e .* (-2i * pi * tt));

29 X2 = sum(b .* e .* (-4 * pi ^ 2 * tt .^ 2));

30 g = 2 * real(conj(X) * X1);

31 h = 2 * real(conj(X1) * X1 + conj(X) * X2);

32 if h >= 0

33 break;

34 end

35 step = -g / h;

36 if abs(step) > df

37 step = sign(step) * df;

38 end

39 f = f + step;

40 if abs(step) < 1e-13 * f

41 break;

42 end

43end

44end