Implementation of ARYule

The graphs show an expt testing  ARYule 100 times on a 2-pole system with 2^8 data points:

>> for i=1:100
AAaryule
drawnow
end
>>

% Aryule method as implemented in AA Book(ref)

N=2^7;n=0:N-1;g=randn(size(n));%3.5*(rand(size(n))-0.5);
%b=1;a=[1 -.9 .3];
b=1;a=conv([1 -.9 .8],[1 .3 .8]);
x=filter(b,a,g);

maxlag=15;
p=4;
[rx,lags]=xcorr(x,x,maxlag,'biased');
rx=rx(maxlag+1:maxlag+1+p);
R=toeplitz(rx);
R1=inv(R);
b2=1/R1(1,1);b0=sqrt(b2);            %find b(0) squared
a1=b2*R1(2:p+1,1);                      %fin AR coefficients a(1) - a(5)

% Obtain a coeffs using 'inner' matrix
RR=R(1:p,1:p);
a2=inv(RR)*-rx(2:p+1)';
[a1(1:p-1) a2(1:p-1)];
[H,w]=freqz(b,a);[H1,w]=freqz(b,[1 a1']);
figure(1)
plot(w,20*log10(abs(H)),w,20*log10(abs(H1)));hold on
axis([0 pi -20 20]);

figure(2)
th=0:.1:2*pi;plot(cos(th),sin(th));hold on
plot(roots([1 a1']),'xr')
%zplane(1,[1 a1']);hold on

Negative Frequency

Search for discussion on negative frequency came up with this page
Negative Frequency

Signal Processing Toolbox - Measuring the Power of Deterministic Periodic Signals Demo

Signal Processing Toolbox 6.12

Measuring the Power of Deterministic Periodic Signals

Richard says:

'This page is an interesting and important in understanding the concept of power measurement / estimation of discrete time signals in the frequency domain.'



This demonstration will focus on power signals, specifically deterministic periodic signals. Although continuous in time, periodic deterministic signals produce discrete power spectrums. For this reason we will use mean-square (power) to measure the signal's power at a specific frequency.
We will provide two examples of how to measure a signal's average power. The examples will use sine waves and assume a load impedance of 1 Ohm.

http://www.mathworks.com/products/signal/demos.html?file=/products/demos/shipping/signal/deterministicsignalpower.html

Use of TFESTIMATE to plot the frequency response of a system from its I/O data

Input and output data with frequency response estimate using tfestimate the implementation of pwelch for cross power spectum. Zooming in on both is interesting:

The m-file:

N=2^16;n=0:N-1;x=randn(size(n));
b=[0 1];a=[1 -1.8 .9];
y=filter(b,a,x);
subplot(1,2,1)
plot(n,y,n,x)
title('Input and output data; N=2^1^6;')
legend('output, y','input, x')
[H,w]=freqz(b,a);
[Pxyh,F]=tfestimate(x,y,hamming(N),0,N,1);
[Pxyr,F]=tfestimate(x,y,rectwin(N),0,N,1);
subplot(1,2,2)
plot(F,abs(Pxyr),F,abs(Pxyh),w/2/pi,abs(H))
title('Frequency response and estimates using tfestimate')
legend('tfe:Rectwin','tfe:Hamming','Freq Resp of filter')

Non Parametric PSD continued. Correction from previous

The power spectrum now matches the frequency response of the filter.

plot(f,10*log10(Py/2),'g','LineWidth',3)

I divided Py by 2 before plotting.

Demonstration of Non-parametric Power Spectrum Estimation

The graph shows spectrum of data passed through a filter calculated using pwelch MATLAB command to evaluate the periodogram by welches method.


%Design nice filter
Wp = [.2 .7], Ws = [.1 .8];Rp=1,Rs=30;

[N,Wn]=ellipord(Wp,Ws,Rp,Rs);

[b,a]=ellip(N,Rp,Rs,Wn);

N=2^12;n=0:N-1;

% generate and filter data

%x=rand(size(n))*2-1;

x=randn(size(n));

y=filter(b,a,x);

%Periodogram

[Px,f] = pwelch(x,rectwin(N),0,N,1);

[Py,f] = pwelch(y,rectwin(N),0,N,1);

plot(f,10*log10(Px),f,10*log10(Py))

hold on

title('Periodogram rectangular window');xlabel('frequency cycles per sample');ylabel('P')

% Bartlett Rectangular window

[Px,f] = pwelch(x,rectwin(N/16),0,N/16,1);

[Py,f] = pwelch(y,rectwin(N/16),0,N/16,1);

plot(f,10*log10(Py),'g','LineWidth',3)

% Bartlett Hamming window

[Px,f] = pwelch(x,hamming(N/16),0,N/16,1);

[Py,f] = pwelch(y,hamming(N/16),0,N/16,1);

plot(f,10*log10(Py),'c','LineWidth',2)

legend('x:periodogram','y:periodogram','y:Bartlett,rect','y:Bartlett,hamming')

xlabel('frequency cycles per sample');ylabel('P')

[H,w]=freqz(b,a);plot(w/2/pi,20*log10(abs(H)),'r')

axis([0 .5 -80 20]);grid

Correlation using Humusoft

Finally got down to using the new analog boards in room 14. I am focussing on block oriented signal processing using the buffer and frame in/out blocks. They are working quite happily at 10kHz sampling and work well at 20 kHz. I have both i/p and o/p going but still have to work out the exact buffering organisation of the boards. I am readyb to get the cross correlation functions going again. I'll post the code tomorrow.

Cross correlation continued....

I have included a better description of the figure.
The xcorr function is used 'unbiased' as I am simply showing that the cross correlation rxw is not equal to rwx.


% m-file showing that rxw /= (not equal to) rwx
% rxw(m)=0 for m<0,>= 0;
% rwx(m)=0 for m>0, =h(-m) m<= 0;
% This is important when deriving the Yule-Walker equations for an AR
% process
N=2^16;n=0:N-1;w=randn(size(n));
b=[0 1];a=[1 -1.8 .9];
x=filter(b,a,w);
M=150;m=-M:M;rxw=xcorr(x,w,M);rwx=xcorr(w,x,M);
subplot(2,1,1);plot(m,rxw,'.');
title('Crosscorrelation rxw');xlabel('m');ylabel('rxw');grid
subplot(2,1,2);plot(m,rwx,'.');xlabel('m');ylabel('rwx');grid

Trying to do Yule walker from Proakis&Manolakis

% m-file showing that rxw /= rwx

% rxw(m)=0 for m<0,>= 0;

% rwx(m)=0 for m>0, =h(-m) m<= 0;

% This is important when deriving the Yule-Walker equations for an AR

% process

N=2^10;n=0:N-1;w=randn(size(n));

b=1;a=[1 -1 .8];

x=filter(b,a,w);

M=40;m=-M:M;rxw=xcorr(x,w,M);rwx=xcorr(w,x,M);

subplot(2,1,1);plot(m,rxw);subplot(2,1,2);plot(m,rwx)