You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
63 lines
1.9 KiB
63 lines
1.9 KiB
6 years ago
|
% mag_to_phase.m
|
||
|
|
%
|
||
|
|
% David Rowe Sep 2015
|
||
|
|
%
|
||
|
|
% Slighly 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
|
||
|
|
|