Homec4science
Diffusion EMPoWER (master) Browse

R1252/simulation-code/cvx/examples/antenna_array_design/html/spectral_fact.html
de096106547amaster

Last Change
Raw File

spectral_fact.html
View Options

<!DOCTYPE HTML>
<html>
<head>
<meta charset="UTF-8">
<title>Spectral factorization using Kolmogorov 1939 approach</title>
<link rel="canonical" href="http://cvxr.com/cvx/examples/antenna_array_design/html/spectral_fact.html">
<link rel="stylesheet" href="../../examples.css" type="text/css">
</head>
<body>
<div id="header">
<h1>Spectral factorization using Kolmogorov 1939 approach</h1>
Jump to:&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#source">Source code</a>&nbsp;&nbsp;&nbsp;&nbsp;
Text output
&nbsp;&nbsp;&nbsp;&nbsp;
Plots
&nbsp;&nbsp;&nbsp;&nbsp;<a href="../../index.html">Library index</a>
</div>
<div id="content">
<a id="source"></a>
<pre class="codeinput">
<span class="comment">% (code follows pp. 232-233, Signal Analysis, by A. Papoulis)</span>
<span class="comment">%</span>
<span class="comment">% Computes the minimum-phase impulse response which satisfies</span>
<span class="comment">% given auto-correlation.</span>
<span class="comment">%</span>
<span class="comment">% Input:</span>
<span class="comment">% r: top-half of the auto-correlation coefficients</span>
<span class="comment">% starts from 0th element to end of the auto-corelation</span>
<span class="comment">% should be passed in as a column vector</span>
<span class="comment">% Output</span>
<span class="comment">% h: impulse response that gives the desired auto-correlation</span>
<span class="keyword">function</span> h = spectral_fact(r)
<span class="comment">% length of the impulse response sequence</span>
nr = length(r);
n = (nr+1)/2;
<span class="comment">% over-sampling factor</span>
mult_factor = 30; <span class="comment">% should have mult_factor*(n) &gt;&gt; n</span>
m = mult_factor*n;
<span class="comment">% computation method:</span>
<span class="comment">% H(exp(jTw)) = alpha(w) + j*phi(w)</span>
<span class="comment">% where alpha(w) = 1/2*ln(R(w)) and phi(w) = Hilbert_trans(alpha(w))</span>
<span class="comment">% compute 1/2*ln(R(w))</span>
w = 2*pi*[0:m-1]/m;
R = exp( -j*kron(w',[-(n-1):n-1]) )*r;
R = abs(real(R)); <span class="comment">% remove numerical noise from the imaginary part</span>
figure; plot(20*log10(R));
alpha = 1/2*log(R);
<span class="comment">% find the Hilbert transform</span>
alphatmp = fft(alpha);
alphatmp(floor(m/2)+1:m) = -alphatmp(floor(m/2)+1:m);
alphatmp(1) = 0;
alphatmp(floor(m/2)+1) = 0;
phi = real(ifft(j*alphatmp));
<span class="comment">% now retrieve the original sampling</span>
index = find(rem([0:m-1],mult_factor)==0);
alpha1 = alpha(index);
phi1 = phi(index);
<span class="comment">% compute the impulse response (inverse Fourier transform)</span>
h = ifft(exp(alpha1+j*phi1),n);
</pre>
</div>
</body>
</html>
c4science · Help