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