📄 the tfarma toolbox v1_0 tutorial.htm
字号:
TFAR parameters completely describe the second-order properties of the given
data!</SPAN> <BR></P>
<H3>2.2 TFAR Modeling in MATLAB<SPAN style="FONT-WEIGHT: bold"></SPAN></H3>
<P>As an example, we generate a TFAR(3, 2) signal <SPAN
style="FONT-STYLE: italic">x</SPAN>[<SPAN style="FONT-STYLE: italic">n</SPAN>]
with <SPAN style="FONT-STYLE: italic">N</SPAN>= 256 as shown above. TF
parameters <SPAN style="FONT-STYLE: italic">c</SPAN><SMALL><SPAN
style="FONT-STYLE: italic">m, l</SPAN></SMALL> (<SPAN
style="FONT-STYLE: italic">c</SPAN> stands for either <SPAN
style="FONT-STYLE: italic">a</SPAN> or <SPAN
style="FONT-STYLE: italic">b</SPAN>) of order (<SPAN
style="FONT-STYLE: italic">M, L</SPAN>) are stored in a (2<SPAN
style="FONT-STYLE: italic">L</SPAN>+1)x(<SPAN
style="FONT-STYLE: italic">M</SPAN>+1) matrix <BIG><BIG><SPAN
style="FONT-FAMILY: monospace">Cml</SPAN></BIG></BIG>, <BIG><BIG><SPAN
style="FONT-FAMILY: monospace">size(Cml)== [2*L+1, M+1]</SPAN></BIG></BIG>, with
the frequency-lag <SPAN style="FONT-STYLE: italic">l</SPAN> as the row index
(from -<SPAN style="FONT-STYLE: italic">L</SPAN> to <SPAN
style="FONT-STYLE: italic">L</SPAN>) and the time-lag <SPAN
style="FONT-STYLE: italic">m</SPAN> as the column index (from 0 to <SPAN
style="FONT-STYLE: italic">M</SPAN>). The model parameters are<BR></P><SPAN
style="FONT-FAMILY: monospace"></SPAN><PRE>N= 256;<BR><BR>Aml= [<BR> 0 0.2774 + 0.1809i -0.0145 + 0.0404i 0.0541 + 0.0448i<BR> 0 0.1762 + 0.1503i 0.0565 - 0.0074i 0.0807 - 0.0163i<BR> 1.0000 -0.2140 - 0.1005i 0.2221 + 0.1832i -0.1783 - 0.0609i<BR> 0 0.0744 + 0.1909i 0.0529 + 0.1038i 0.1087 - 0.0698i<BR> 0 0.1601 + 0.1545i -0.0889 - 0.0562i -0.0116 + 0.0304i<BR>];<BR><BR>B0l= [<BR> -0.0455 + 0.0222i<BR> 0.0008 - 0.0743i<BR> 0.1617<BR> 0.0008 + 0.0743i<BR> -0.0455 - 0.0222i<BR>];<BR></PRE>
<P>[Plot the parameter functions. ] Note that for <SPAN
style="FONT-STYLE: italic">L</SPAN>= 0 the TFAR(<SPAN
style="FONT-STYLE: italic">M, L</SPAN>) model would collapse to a usual AR(<SPAN
style="FONT-STYLE: italic">M</SPAN>) model with parameters </P><PRE><SPAN style="FONT-FAMILY: monospace">Aml= [</SPAN>1.0000 -0.2140 - 0.1005i 0.2221 + 0.1832i -0.1783 - 0.0609i<SPAN style="FONT-FAMILY: monospace">];<BR><BR></SPAN><SPAN style="FONT-FAMILY: monospace">B0l= [</SPAN>0.1617<SPAN style="FONT-FAMILY: monospace">];</SPAN><BR></PRE>
<P>A TFAR model is thought to generate the nonstationary signal <SPAN
style="FONT-STYLE: italic">x</SPAN>[<SPAN style="FONT-STYLE: italic">n</SPAN>]
by passing normalized stationary white Gaussian noise <SPAN
style="FONT-STYLE: italic">e</SPAN>[<SPAN style="FONT-STYLE: italic">n</SPAN>]
(<BIG><BIG><SPAN style="FONT-FAMILY: monospace">randn(N, 1)</SPAN></BIG></BIG>)
through a time-varying IIR filter. The filter parameters <SPAN
style="FONT-STYLE: italic">c<SMALL>m</SMALL></SPAN>[<SPAN
style="FONT-STYLE: italic">n</SPAN>] are generated by a width-<SPAN
style="FONT-STYLE: italic">L</SPAN> Fourier series with the model parameters
<SPAN style="FONT-STYLE: italic">c<SMALL>m, l</SMALL></SPAN>. A process
realization can be generated by (<BIG><BIG><SPAN
style="FONT-FAMILY: monospace">help </SPAN><SPAN
style="FONT-FAMILY: monospace">tfarma_gen</SPAN><SPAN
style="FONT-FAMILY: monospace"></SPAN></BIG></BIG>) <BR></P><PRE><SPAN style="FONT-FAMILY: monospace">x= tfarma_gen(randn(N, 1), Aml, B0l);</SPAN></PRE>
<P><SPAN style="FONT-STYLE: italic"></SPAN>The evolutionary TFAR(<SPAN
style="FONT-STYLE: italic">M</SPAN>) spectrum (expected quantity) can be
computed from the given TF parameters <BIG><BIG><SPAN
style="FONT-FAMILY: monospace">Aml</SPAN><SMALL><SMALL>, </SMALL></SMALL><SPAN
style="FONT-FAMILY: monospace">B0l</SPAN></BIG></BIG> as <BR></P><PRE><SPAN style="FONT-FAMILY: monospace">Px= tfarma_evsp(Aml, B0l, N);</SPAN></PRE>
<P>It looks like (compare this spectrum to (<BIG><BIG><SPAN
style="FONT-FAMILY: monospace">figure(2)</SPAN></BIG></BIG>!) <BR></P><PRE><SPAN style="FONT-FAMILY: monospace">figure(4);clf;imagesc(rot90(log(Px)));colormap(flipud(gray));<BR></SPAN><SPAN style="FONT-FAMILY: monospace">set(gca, 'XTick', [1 N/4 N/2 3*N/4 N])<BR>set(gca, 'XTickLabel', [0 N/4-1 N/2-1 3*N/4-1 N-1])<BR>set(gca, 'YTick', [1 N/4 N/2 3*N/4 N])<BR>set(gca, 'YTickLabel', [N/2-1 N/4-1 0 -N/4 -N/2])<BR><BR></SPAN></PRE>
<P></P><PRE><SPAN style="FONT-FAMILY: monospace"></SPAN></PRE><IMG
style="WIDTH: 75%; HEIGHT: 50%" alt="Expected evolutionary TFAR spectrum"
src="The TFARMA Toolbox v1_0 Tutorial_files/sigspec.jpg"><BR>
<H3>2.3 Estimator Accuracy</H3>
<P>The estimated model order (<BIG><BIG><SPAN
style="FONT-FAMILY: monospace">MEst</SPAN></BIG></BIG>, <BIG><BIG><SPAN
style="FONT-FAMILY: monospace">LEst</SPAN></BIG></BIG>) need not to be the same
as the given one (3, 2). Repeat the following lines until a model of order (3,
2) is found! <BR></P><PRE><SPAN style="FONT-FAMILY: monospace"></SPAN><SPAN style="FONT-FAMILY: monospace">[AmlEst, B0lEst]= tfar(</SPAN><SPAN style="FONT-FAMILY: monospace">tfarma_gen(randn(N, 1), Aml, Bml));<BR></SPAN><SPAN style="FONT-FAMILY: monospace">[MEst, LEst]= param_dim(</SPAN><SPAN style="FONT-FAMILY: monospace">AmlEst</SPAN><SPAN style="FONT-FAMILY: monospace">)</SPAN><BR></PRE>
<P>Now, we are able to compute how much the estimated parameters deviate from
the given ones. We simply compute the norm of the difference <BR></P><PRE>norm([Aml B0l] - [AmlEst B0lEst])<BR></PRE>
<P>Further, the estimated model order should converge to the true order. If you
have some time left, run <BR></P><PRE>ORDER= [];<BR>for mm= 1:100<BR> mm<BR> [AmlEst, B0lEst]= tfar(tfarma_gen(randn(N, 1), Aml, Bml));<BR> <SPAN style="FONT-FAMILY: monospace">[MEst, LEst]= </SPAN>param_dim(AmlEst);<BR> ORDER= [ORDER; [MEst LEst]];<BR>end;<BR>mean(ORDER)<BR></PRE>
<H3>2.4 <SPAN style="FONT-WEIGHT: bold"></SPAN>Example-MATLAB-Script<SPAN
style="FONT-WEIGHT: bold"></SPAN></H3>
<P>The MATLAB script given below performs classical (nonparametric) TF-spectral
analysis of the given [only for complex- valued] nonstationary data <SPAN
style="FONT-STYLE: italic">x</SPAN>[<SPAN style="FONT-STYLE: italic">n</SPAN>]
and it also fits a TFAR model. First, the mean is subtracted and two histograms
(real and imaginary) are plotted into <BIG><BIG><SPAN
style="FONT-FAMILY: monospace">figure(1)</SPAN></BIG></BIG>. To embed the given
signal of length <SPAN style="FONT-STYLE: italic">N</SPAN><FONT size=-1>0</FONT>
into a power-of-two frame of length <SPAN style="FONT-STYLE: italic">N</SPAN>,
we generate blocks of Gaussian noise (SNR= 10dB) which are merged before and
after the given signal. Then, the modified signal <SPAN
style="FONT-STYLE: italic">y</SPAN>[<SPAN style="FONT-STYLE: italic">n</SPAN>]
(<BIG><BIG><SPAN style="FONT-FAMILY: monospace">figure(2)</SPAN></BIG></BIG>),
its magnitude squared spectrogram, i.e., physical spectrum (<BIG><BIG><SPAN
style="FONT-FAMILY: monospace">figure(3)</SPAN></BIG></BIG>), the Rihaczek
distribution with the same windowing as for the spectrogram (<BIG><BIG><SPAN
style="FONT-FAMILY: monospace">figure(4)</SPAN></BIG></BIG>) and the estimated
evolutionary TFAR spectrum (<BIG><BIG><SPAN
style="FONT-FAMILY: monospace">figure(5)</SPAN></BIG></BIG>) are shown. Note
that <SPAN style="FONT-STYLE: italic">N</SPAN> is limited to ca. 1024-2048
(depends on your computer power). <BR></P><PRE>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<BR>% the variable x needs to contain your complex-valued <BR>% signal of any length<BR>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<BR>% Into a column vector<BR>x= x(:);<BR><BR>% Mean subtraction<BR>x= x-mean(x);<BR><BR>% Real block length<BR>N0= length(x);<BR><BR>% Histogram plot<BR>figure(1);clf;<BR>subplot(1, 2, 1);hist(real(x));colormap(gray);<BR>subplot(1, 2, 2);hist(imag(x));colormap(gray);<BR><BR>% Embedding in a power-of-two-frame with SNR= 10dB<BR>n= ceil(log2(N0));<BR>N= 2^n;<BR>varre= var(real(x));<BR>varim= var(imag(x));<BR><BR>v1= sqrt(varre/10)*randn(ceil ((N-N0)/2), 1) + j*sqrt(varim/10)*randn(ceil ((N-N0)/2), 1);<BR>v2= sqrt(varre/10)*randn(floor((N-N0)/2), 1) + j*sqrt(varim/10)*randn(floor((N-N0)/2), 1);<BR><BR>y= [v1; x; v2];<BR>maxx= .1*ceil(max(abs(y))*11)<BR><BR>% Embedded signal plot<BR>figure(2);clf;<BR>subplot(2, 1, 1);plot(real(y), 'Color', 'k', 'Linewidth', 2);axis([1 N -maxx maxx]);<BR>set(gca, 'XTick', [1 N/4 N/2 3*N/4 N])<BR>set(gca, 'XTickLabel', [0 N/4-1 N/2-1 3*N/4-1 N-1])<BR>subplot(2, 1, 2);plot(imag(y), 'Color', 'k', 'Linewidth', 2);axis([1 N -maxx maxx]);<BR>set(gca, 'XTick', [1 N/4 N/2 3*N/4 N])<BR>set(gca, 'XTickLabel', [0 N/4-1 N/2-1 3*N/4-1 N-1])<BR><BR>% The FFT length for the spectrogram<BR>NFFT= N/4;<BR><BR>% The window<BR>w= symmpad(hanning(2*ceil(NFFT/6)-1), NFFT);<BR><BR>% Spectrogram<BR>S= specgram([ zeros(N/8, 1); y; zeros(N/8-1, 1)], NFFT, [], w, NFFT-1);<BR>S= [S(N/8+1:end, :); S(1:N/8, :)];<BR>figure(3);imagesc(log(flipud(abs(S).^2)));colormap(flipud(gray));<BR>set(gca, 'XTick', [1 N/4 N/2 3*N/4 N])<BR>set(gca, 'XTickLabel', [0 N/4-1 N/2-1 3*N/4-1 N-1])<BR>set(gca, 'YTick', [1 N/16 N/8 3*N/16 N/4])<BR>set(gca, 'YTickLabel', [N/2-1 N/4-1 0 -N/4 -N/2])<BR>drawnow<BR><BR>% Rihaczek distribution<BR>Psi= tf_multiwin(N, 16, 16);<BR>WxEst= real(ml_to_nk( ambi_est(y, y, -1).*conj(Psi) ));<BR>figure(4);clf;imagesc(rot90(log(WxEst)));colormap(flipud(gray));<BR>set(gca, 'XTick', [1 N/4 N/2 3*N/4 N])<BR>set(gca, 'XTickLabel', [0 N/4-1 N/2-1 3*N/4-1 N-1])<BR>set(gca, 'YTick', [1 N/4 N/2 3*N/4 N])<BR>set(gca, 'YTickLabel', [N/2-1 N/4-1 0 -N/4 -N/2])<BR><BR>% TFAR Fit<BR>[AmlEst, BmlEst]= tfar(y);<BR>[MEst, LEst]= param_dim(AmlEst)<BR>PxEst= tfarma_evsp(AmlEst, BmlEst, N);<BR>figure(5);clf;imagesc(log(rot90(PxEst, 1)));colormap(flipud(gray));<BR>set(gca, 'XTick', [1 N/4 N/2 3*N/4 N])<BR>set(gca, 'XTickLabel', [0 N/4-1 N/2-1 3*N/4-1 N-1])<BR>set(gca, 'YTick', [1 N/4 N/2 3*N/4 N])<BR>set(gca, 'YTickLabel', [N/2-1 N/4-1 0 -N/4 -N/2])<BR></PRE>
<P></P>
<P><SPAN style="FONT-WEIGHT: bold">At this point, you should be able to apply
the TFARMA Toolbox in its simplest way (just type </SPAN><BIG
style="FONT-WEIGHT: bold"><BIG><SPAN
style="FONT-FAMILY: monospace"><BR></SPAN></BIG></BIG></P><PRE><BIG style="FONT-WEIGHT: bold"><BIG><SPAN style="FONT-FAMILY: monospace">[AmlEst, B0lEst]= tfar(x)</SPAN></BIG></BIG><SPAN style="FONT-WEIGHT: bold"> </SPAN></PRE><PRE><SPAN style="FONT-WEIGHT: bold"></SPAN></PRE>
<P><SPAN style="FONT-WEIGHT: bold">with <SPAN
style="FONT-STYLE: italic">N</SPAN> a power of two</SPAN><SPAN
style="FONT-WEIGHT: bold">)</SPAN>. You can now fit a TFAR(<SPAN
style="FONT-STYLE: italic">M, L</SPAN>) model to your nonstationary data. The
estimated model parameters (<SPAN
style="FONT-STYLE: italic">a</SPAN><SMALL><SPAN style="FONT-STYLE: italic">m,
l</SPAN></SMALL><SPAN style="FONT-STYLE: italic">, b</SPAN><SMALL>0<SPAN
style="FONT-STYLE: italic">, l</SPAN></SMALL>) serve as a complete statistical
second-order description of the given nonstationary process <SPAN
style="FONT-STYLE: italic">x</SPAN>[<SPAN style="FONT-STYLE: italic">n</SPAN>]
and further processing of <SPAN style="FONT-STYLE: italic">x</SPAN>[<SPAN
style="FONT-STYLE: italic">n</SPAN>] can be done based on the estimated
reduced-detail parameter set<SPAN style="FONT-STYLE: italic"></SPAN><SMALL><SPAN
style="FONT-STYLE: italic"></SPAN></SMALL><SPAN
style="FONT-STYLE: italic"></SPAN><SMALL><SPAN
style="FONT-STYLE: italic"></SPAN></SMALL>. <BR><SPAN
style="FONT-WEIGHT: bold"></SPAN></P>
<P><SPAN style="FONT-WEIGHT: bold"></SPAN><SPAN
style="FONT-WEIGHT: bold"></SPAN></P>
<H3>2.5 Vector TFAR Models<BR></H3>
<P>If you are given a multivariate (vector) signal <SPAN
style="FONT-WEIGHT: bold; FONT-STYLE: italic">x</SPAN>[<SPAN
style="FONT-STYLE: italic">n</SPAN>] with <SPAN
style="FONT-STYLE: italic">F</SPAN> components <SPAN
style="FONT-STYLE: italic">x<SMALL>1</SMALL></SPAN>[<SPAN
style="FONT-STYLE: italic">n</SPAN>] to <SPAN
style="FONT-STYLE: italic">x<SMALL>F</SMALL></SPAN>[<SPAN
style="FONT-STYLE: italic">n</SPAN>], you will enjoy fitting a VTFAR model. The
MATLAB syntax stays the same, just write a <BIG><BIG><SPAN
style="FONT-FAMILY: monospace">V</SPAN></BIG></BIG> in front of
<BIG><BIG><SPAN style="FONT-FAMILY: monospace">TFARMA</SPAN></BIG></BIG>! The
model parameter matrices <SPAN
style="FONT-WEIGHT: bold; FONT-STYLE: italic">c</SPAN><SMALL
style="FONT-STYLE: italic">m, l</SMALL> are stored as 4D arrays <BIG><BIG><SPAN
style="FONT-FAMILY: monospace">CML</SPAN></BIG></BIG>, <BIG><BIG><SPAN
style="FONT-FAMILY: monospace">size(CML)== [F, F, 2*L+1,
M+1]</SPAN></BIG></BIG>. The signal is stored in an <SPAN
style="FONT-STYLE: italic">N</SPAN>x<SPAN style="FONT-STYLE: italic">F</SPAN>
matrix <BIG><BIG><SPAN style="FONT-FAMILY: monospace">X</SPAN></BIG></BIG>
(instead of <BIG><BIG><SPAN
style="FONT-FAMILY: monospace">x</SPAN></BIG></BIG>). Similar to the scalar
case, we simply type (note that we write <BIG><BIG><SPAN
style="FONT-FAMILY: monospace">AMLEst</SPAN></BIG></BIG> instead of
<BIG><BIG><SPAN style="FONT-FAMILY: monospace">AmlEst</SPAN></BIG></BIG>!)
<BR></P><PRE>[AMLEst, B0LEst]= vtfar(X)<BR><SPAN style="FONT-FAMILY: monospace">[MEst, LEst]= param_dim(</SPAN><SPAN style="FONT-FAMILY: monospace">AMLEst</SPAN><SPAN style="FONT-FAMILY: monospace">)<BR></SPAN><SPAN style="FONT-FAMILY: monospace"></SPAN></PRE>
<P>we obtain a [nonstabilized] VTFAR(<SPAN style="FONT-STYLE: italic">M,
L</SPAN>) model. To generate an example signal <SPAN
style="FONT-WEIGHT: bold; FONT-STYLE: italic">x</SPAN>[<SPAN
style="FONT-STYLE: italic">n</SPAN>], we choose a monic VTFAR(2, 1) with <SPAN
style="FONT-STYLE: italic">F</SPAN>= 2, we type <BR></P><PRE>N= 256;<BR>AML(:,:,1,1)= [ 0 0<BR> 0 0];<BR>AML(:,:,2,1)= [ 1 0<BR> 0 1];<BR>AML(:,:,3,1)= [ 0 0<BR> 0 0];<BR>AML(:,:,1,2)= [ -0.0301 + 0.0801i 0.0422 + 0.1571i<BR> -0.0488 - 0.0455i 0.1999 + 0.0740i];<BR>AML(:,:,2,2)= [ 0.2794 - 0.0079i 0.1020 - 0.0698i<BR> -0.1316 - 0.1935i 0.3989 + 0.2268i];<BR>AML(:,:,3,2)= [ -0.0735 - 0.2695i -0.2218 + 0.1556i<BR> 0.2284 - 0.0626i -0.0291 + 0.1997i];<BR>AML(:,:,1,3)= [ -0.1848 + 0.1263i 0.0352 + 0.2071i<BR> -0.0592 + 0.0581i 0.0468 + 0.1225i];<BR>AML(:,:,2,3)= [ 0.1046 + 0.0451i -0.0107 + 0.0301i<BR> 0.0745 + 0.0095i 0.0314 - 0.0086i];<BR>AML(:,:,3,3)= [ -0.0042 + 0.0378i -0.0433 + 0.1037i<BR> -0.0274 - 0.2137i -0.0920 + 0.1374i];<BR>B0L= eye(2);<BR>X= vtfarma_gen(randn(N, 2), AML, B0L);<BR></PRE>
<P></P>
<H2>[3. Second Step: TFARMA Identification]</H2>
<P>[not available yet] For a first try, TFAR identification might be sufficient
for your application. However, if the signal contains not only spectral peaks,
but also deep spectral notches, a TFMA or a TFARMA model would fit the data much
better. <BR></P>
<H3>[3.1 Vector-TFARMA Identification]</H3>
<P>[not available yet] If a multivariate nonstationary signal needs to be
analyzed, one can use the VTFARMA method presented here. <BR></P>
<ADDRESS><SPAN style="FONT-WEIGHT: bold"></SPAN><SPAN
style="FONT-WEIGHT: bold"></SPAN><SPAN
style="FONT-WEIGHT: bold"></SPAN>Questions, comments: <A
href="mailto:michael.jachan@ieee.org">michael.jachan at
ieee.org</A>
last
changes: 301204</ADDRESS>
<ADDRESS><BR></ADDRESS></BODY></HTML>
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -