Appendix 7

correlation.m: MATLAB code for calculating filtered and unfiltered correlation between experimental data and RAYNOISE generated impulse responses.

Permits optional filtering with parameter specification for either Butterworth or Chebyshev type II low-pass filters for correlation analysis

clear

%	CORRELATION.M (NEWANAL.M)
%	calculates correlation between resampled experimental impulses
%	and RAYNOISE generated .wav file impulse responses. Allows options
%	for filtering, including either butterworth or chebyshev type II
%	low-pass filters


% Allow for variable user input

part1='expimp';
part2='pt';
part3='impm';

part5=input('Which model are you analyzing? >>> ','s');
part4=input('Which point are you analyzing? >>> ','s');
part6='.wav';

% determine filenames

expfilename=[part1 part4]
calcfilename=[part2 part4 part3 part5 part6]

%load data from .mat and RAYNOISE generated .wav files

load(expfilename);
[y,fs,format]=wread16(calcfilename);
fs=fs
fnyquist=fs*0.5;
right=y(1,:);
left=y(2,:);

%data=data./max(abs(data));
data=data-mean(data);
samps=length(data);

startsamp=1
endsamp=round(samps*(1))+startsamp-1
samps=endsamp-startsamp;

% data filtering section

filteranswer=input('Do you wish to filter the data? >>> ','s');

filterokay='n';
if filteranswer=='y',
	while filterokay ~= 'y',
		fcutoff=input('Enter cutoff frequency >>> ');
		normcutoff=fcutoff/fnyquist	% butter needs a cuttoff freq
		between 0 and 1
		order=input('Enter a filter order >>> ');
		fftpoints=1024;
		type=input('Enter A for chebychev II or B for Butterworth type
		filter. >>> ','s');
		if type=='a',
			dbripple=input('Enter any acceptable stopband dB ripple
			>>> ');
			[b,a]=cheby2(order,dbripple,normcutoff);
		end
		if type=='b',
			[b,a]=butter(order,normcutoff);
		end
		[h,w]=freqz(b,a,fftpoints);
		figure(1);
		freqz(b,a,fftpoints);title('data filter');
		print -dps filter.ps;
		filterokay=input('is the filter design okay? >>> ','s');
	end
	disp('filter design over');
	disp('filtering data...');
	filtdata=filter(b,a,data);
	left=filter(b,a,left);
	right=filter(b,a,right);
	figure(2);clg;subplot(3,1,1);plot(data);ylabel('dat');title('exp. data')
	subplot(3,1,2);plot(filtdata);ylabel('filt. dat.')
	echogram=abs(filtdata)./max(abs(filtdata));clear filtdata
	subplot(3,1,3);plot(echogram);ylabel('echogram');xlabel('samples @ 44.1
	kHz');
	subplot(3,1,1);title('Filtered experimental data and normalized
	echogram');
	print -dps filter.ps
else 
	filtdata=data;
	figure(2);clg;subplot(2,1,1);plot(data);ylabel('dat');title('exp. data')
	echogram=abs(filtdata)./max(abs(filtdata));clear filtdata
	subplot(2,1,2);plot(echogram);ylabel('echogram');xlabel('samples @ 44.1
	kHz');
	subplot(2,1,1);title('Unfiltered experimental data and normalized
	echogram');
	part9='unfiltered';
end


lengthright=length(right);
lengthleft=length(left);

%normalize RAYNOISE data

normright=abs(right)./max(abs(right));
normleft=abs(left)./max(abs(left));
sum=normright+normleft;
normsum=abs(sum)./max(abs(sum));

figure(3);
if filteranswer=='y',
	subplot(3,1,1);plot(normright);ylabel('right');
	subplot(3,1,2);plot(normleft);ylabel('left');
	subplot(3,1,3);plot(normsum);ylabel('norm sum');	
	xlabel('samples @ 44.1 kHz')
	subplot(3,1,1);title('Normalized filtered RAYNOISE echograms, left,
	right superimposed')
	part9='filtnorm';
	%printfilename=[part2 part4 part3 part5 part9 part8]
	%print -dps [part2 part4 part3 part5 part9 part8]
else
	subplot(3,1,1);plot(normright);ylabel('right');
	subplot(3,1,2);plot(normleft);ylabel('left');
	subplot(3,1,3);plot(normsum);ylabel('norm sum');	
	xlabel('samples @ 44.1 kHz')
	subplot(3,1,1);title('Normalized unfiltered RAYNOISE echograms, left,
	right superimposed')
	part9='filtnorm';
end

disp('calculationg correlations');

% correct for minor length differences between experimental and RAYNOISE data

lengths=[lengthright lengthleft];
echogram=echogram(1:min(lengths));

% calculate correlations and normalize

rightcorr=xcorr(normright,echogram,'coeff');
leftcorr=xcorr(normleft,echogram,'coeff');
sumcorr=xcorr(normsum,echogram,'coeff');

rightcorr=fftshift(rightcorr);
leftcorr=fftshift(leftcorr);
sumcorr=fftshift(sumcorr);

figure(4);clg
subplot(3,1,1);plot(rightcorr(1:samps*0.05));ylabel('rightcorr');
subplot(3,1,2);plot(leftcorr(1:samps*0.05));ylabel('leftcorr');
subplot(3,1,3);plot(sumcorr(1:samps*0.05));ylabel('sumcorr');
xlabel('samples @ 44.1 kHz')
if filteranswer=='y',
	subplot(3,1,1);title('Correlation of filtered experimental and RAYNOISE
	echograms')
	part9='filtcorrs';
	part3='m';
else
	subplot(3,1,1);title('Correlation of unfiltered experimental and
	RAYNOISE echograms')
	part9='unfiltcorrs';
	part3='m';
	
end