1% mag_to_phase.m
2%
3% David Rowe Sep 2015
4%
5% Slighly modified version of http://www.dsprelated.com/showcode/20.php
6%
7% Given a magnitude spectrum in dB, returns a minimum-phase phase
8% spectra.  Both must be sampled at a Nfft. My understanding of this
9% is rather dim, but a working example is good place to start!
10
11
12function [phase s] = mag_to_phase(Gdbfk, Nfft = 512, verbose_en = 0)
13
14  Ns = length(Gdbfk); if Ns~=Nfft/2+1, error("confusion"); end
15  Sdb = [Gdbfk,Gdbfk(Ns-1:-1:2)]; % install negative-frequencies
16
17  S = 10 .^ (Sdb/20); % convert to linear magnitude
18  s = ifft(S);        % desired impulse response
19  s = real(s);        % any imaginary part is quantization noise
20  tlerr = 100*norm(s(round(0.9*Ns:1.1*Ns)))/norm(s);
21  if verbose_en
22    disp(sprintf(['  Time-limitedness check: Outer 20%% of impulse ' ...
23               'response is %0.2f %% of total rms'],tlerr));
24  end
25  % = 0.02 percent
26
27  if verbose_en
28    if tlerr>1.0 % arbitrarily set 1% as the upper limit allowed
29      disp('  Increase Nfft and/or smooth Sdb\n');
30    end
31  end
32
33  c = ifft(Sdb); % compute real cepstrum from log magnitude spectrum
34
35  % Check aliasing of cepstrum (in theory there is always some):
36
37  caliaserr = 100*norm(c(round(Ns*0.9:Ns*1.1)))/norm(c);
38  if verbose_en
39    disp(sprintf(['  Cepstral time-aliasing check: Outer 20%% of ' ...
40                 'cepstrum holds %0.2f %% of total rms\n'],caliaserr));
41  end
42
43  if verbose_en
44    if caliaserr>1.0 % arbitrary limit
45      disp('  Increase Nfft and/or smooth Sdb to shorten cepstrum\n');
46    end
47  end
48
49  % Fold cepstrum to reflect non-min-phase zeros inside unit circle:
50
51  cf = [c(1), c(2:Ns-1)+c(Nfft:-1:Ns+1), c(Ns), zeros(1,Nfft-Ns)];
52
53  Cf = fft(cf); % = dB_magnitude + j * minimum_phase
54
55  % The maths says we are meant to be using log(x), not 20*log10(x),
56  % so we need to scale the phase to account for this:
57  % log(x) = 20*log10(x)/scale;
58
59  scale = (20/log(10));
60  phase = imag(Cf)/scale;
61endfunction
62
63