Measurement of speed of sound OR measurement of distance using sound

This is set up as a lab experiment demonstrating the use of cross correlation. The PC sound card is used together with wavplay  and  wavrecord in MATLAB to generate and measure, using two microphones a fixed distance apart, a noise burst. Wavplay is used in asynchronous mode to play the noise an immediately return control to the command window so that the sound can be recoded using  wavrecord.

Fs=44100;
T=1/Fs;
s=.1;
N=s*Fs;
x=idinput(N);%,'prbs');
xx=zeros(N,2);xx(:,2)=x;
t=0:1/Fs:(N-1)/Fs;
wavplay(xx,Fs,'async')
y=wavrecord(s*Fs,Fs,2);
figure(1)
subplot(2,1,1)
plot(t,y);%plot(t,y(:,1))
subplot(2,1,2)
plot(t,y(:,2))
figure(2)
%subplot(2,1,2)
rxx=xcorr(y(:,1),y(:,2),200);
l=-200:200;
plot(l,rxx,'r');grid on
[A,I]=max(rxx);
Tlag=T*(200-I);
v=340;
d=v*Tlag   % Display distance in command window

Demonstration of Ensemble Stochastic Process

 % Demonstrate ensemble stochastic process
M=4000;N=50;Nsub=6;
x=randn(M,N);

% Uncomment next two lines when we filter the process
% b=1;a=[1 -1.5 .85];
% x=filter(b,a,x')';

n=0:N-1;
nn=20;mm=21;   %  nn and mm are the positions of the markers
for i=1:Nsub
subplot(Nsub,1,i);
stem(n,x(i,:),'MarkerSize',3);%axis([0 N -2.5 2.5]);
hold on
plot(nn,x(i,nn+1),'o','MarkerSize',6,'MarkerFaceColor','r','MarkerEdgeColor','k')
plot(mm,x(i,mm+1),'^','MarkerSize',6,'MarkerFaceColor','g','MarkerEdgeColor','k')
hold off
end

ac=mean(x(:,nn).*x(:,mm))
subplot(Nsub,1,1)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

title(['r_{xx}(',num2str(nn),',',num2str(mm),')  =  ',num2str(ac)],'FontSize',15)

Visualising the Complex Exponential

 

x=-20:.1:20;y=-20:.1:20;
[X,Y] = meshgrid(x,y);
a=-0.05;b=0.75;
% Real Component :    e^(ax)cos(by)
z=exp(a*X).*cos(b*Y);
subplot(2,2,1)
mesh(X,Y,z)
zrr=exp(a*x).*cos(b*0);
hold on
plot3(x,zeros(length(x)),zrr,'r','LineWidth',3)
zir=exp(a*0).*cos(b*y);%zi=exp(j*b*y);
plot3(zeros(length(x)),y,zir,'r','LineWidth',4)
zrr=exp(a*x).*cos(b*y);%zrr.*zir;
plot3(x,y,zrr,'b','LineWidth',4)
hold off

% Imaginary Component :    je^(ax)sin(by)
z=exp(a*X).*sin(b*Y);
subplot(2,2,2)
mesh(X,Y,z)
hold on
zri=exp(a*x)*sin(0);
plot3(x,zeros(length(x)),zri,'r','LineWidth',3)
zii=exp(a*0)*sin(b*y);
plot3(zeros(length(x)),y,zii,'r','LineWidth',4)
zii=exp(a*x).*sin(b*y);
plot3(x,y,zii,'b','LineWidth',4)
hold off

subplot(2,2,3)
plot(zrr+j*zii)
 

3-D view of a waveform

I have used the DFT to decompose the waveform into its components and have used 'plot3' to display the in a 3-D image. The image can be rotated in MATLAB to examine the frequency domain as well as the time domain. The second figure shows the figure rotated.

N=128;n=0:N-1;fund=3/N;%0+.0167/2;
x=sin(2*pi*fund*n)+sin(2*3*pi*fund*n)/3+...
sin(2*pi*5*fund*n)/5+sin(2*pi*7*fund*n)/7;
%x=x.*(hamming(N))';
DFT3D_1(x)

function DFT3D_1(x)
N=length(x);
n=0:N-1;
%x=x.*(hamming(N))';
X=fft(x)*2/N;
X1=abs(X')*ones(1,N);X1=X1';

f=(0:N-1)/N; % Harmonics
[ff,nn]=meshgrid(f,n);
W=cos(2*pi*f'*n);

X1=W'.*X1;
X1(:,1)=x;

plot3(ff,nn,X1,'LineWidth',1);grid on

xlabel('frequency');
ylabel('time seconds');
axis([0 1 0 N -1.3 1.3]);
    view([37.5 30])

Yule Walker Method Versus Welch Method

In this section I compare ARYule with Pwelch for a short data set. The emphasis is on the superior performance of the Yule Method for the short data set.

The figure shows a comparison between the YULE method of PSD Estimation and Welches method for the same amount of data. The code is shown below. Comments please. I can write a lot about this but it is self explanatory. DO you think?

ARYule Figure:
b=1;a=[1 -0.9 .5];
[H,w]=freqz(b,a);

for i=1:1
N=2^7;n=0:N-1;w=randn(size(n));

x=filter(b,a,w);

%ARYuleRGH1;
A=aryule(x,2);

[H1,w]=freqz(b0,A);
subplot(1,2,1)
plot(w,10*log10(abs(H)))
hold on
plot(w,10*log10(abs(H1)),'r')
subplot(1,2,2)
th=0:.05:2*pi;plot(cos(th),sin(th));hold on;
plot(roots(a),'x')
plot(roots(A),'xr')
drawnow
end

subplot(1,2,1);plot(w,10*log10(abs(H)),'LineWidth',3)
title('Yule-Walker Method: N=2^7,M=2^6, no overlap')
hold off
subplot(1,2,2);plot(roots(a),'x','LineWidth',3)
title('POLE plot')
hold off




PWELCH Figure

b=1;a=[1 -0.9 .5];
[H,w]=freqz(b,a);
for i=1:1
N=2^7;n=0:N-1;v=randn(size(n));M=2^6;
x=filter(b,a,v);
[Px,F]=pwelch(x,rectwin(M),0,M,1);
plot(F,10*log10(Px/2),'r')
hold on
drawnow
end
plot(w/2/pi,10*log10(abs(H).^2),'LineWidth',3)
title('welch method: N=2^7,M=2^6, no overlap')
hold off

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