Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- % Modified by: Brian Cai
- % ME 622: Advanced Acoustics
- % Assignment 2: Problem 2
- % Due September 18, 2014
- % This m file computes the autospectrum, prob density, and transfer
- % funtions from a time domain signal acquired by the computer. It is
- % written using Matlab Version 7.0. Note that to acquire data using this
- % progam the user must make sure that the desired recording device is
- % "selected" from the "Record Control" dialog box in Windows. To do this,
- % right click on the loudspeaker icon in the quick launch section at the
- % right end of the task bar and select "Open Volume Control." Click
- % "options" and select "properties" and then select "recording". Then
- % select "Line-In" or "Microphone" or "CD Audio" etc. depending on which
- % input you want to record from. If you don't have the loudspeaker icon in
- % your quick launch section you can select the recording devices from
- % Control Panel>Sounds and Audio Devices.
- % R. N. Miles
- clear all; close all; clc
- % The number of records to acquire.
- numrecord = 50;
- % num is the total number of time points per record. This should be a power
- % of 2 for the FFT algorithm.
- num = 16384;
- % The sample rate in samples/second. This is a reasonable rate for audio
- % signals.
- fs = 44100;
- % fs = 96000;
- % The amount of time in each record.
- Trecord = num / fs;
- % Total time in seconds.
- totaltime = numrecord * Trecord;
- % Acquire data from the sound card at the rate fs (samples/s) using 16
- % bits, two channels.
- recorder = audiorecorder(fs , 16 , 2);
- % recorder = audiorecorder(fs , 24 , 2);
- recordblocking(recorder , totaltime);
- audioarray = getaudiodata(recorder);
- % The rest of the code remains unchanged.
- delfreq=1/Trecord; %The frequency spacing in Hertz.
- numt=length(audioarray(:,1));
- time=linspace(0,Trecord,num); %Time for first record
- figure ('name','time domain plot')
- plot(time,audioarray(1:num,1),time,audioarray(1:num,2))
- xlabel('Time in seconds')
- ylabel('original signal')
- title('Record 1 of original signal')
- grid on
- delt=1/fs;
- freq=linspace(0,fs,num);
- spec1=0*freq'; %Create an array of zeros that is the same size as freq.
- spec2=spec1;
- crossspec=spec1; %Create an array of zeros that is the same size as freq.
- ex2time=0;
- for irec=1:numrecord
- %rtime1=audioarray(1+num*(numrecord-1):num+num*(numrecord-1),1); %Store the first num points from channel 1 in an array
- %rtime2=audioarray(1+num*(numrecord-1):num+num*(numrecord-1),2); %Store the first num points from channel 2 in an array
- rtime1=audioarray(1+num*(irec-1):num+num*(irec-1),1); %Store the first num points from channel 1 in an array
- rtime2=audioarray(1+num*(irec-1):num+num*(irec-1),2); %Store the first num points from channel 2 in an array
- fft1=fft(rtime1,num);
- fft2=fft(rtime2,num);
- ex2time=ex2time+sum(rtime1.^2)/num; %Mean square of each time record
- spec1=spec1+abs(fft1).^2*2*delt/num; %calculate the autospectrum of the random signal.
- spec2=spec2+abs(fft2).^2*2*delt/num; %calculate the autospectrum of the random signal.
- % This will be a single-sided spectrum in units^2/Hz
- crossspec=crossspec+fft1.*conj(fft2)*2*delt/num;
- end
- spec1=spec1/numrecord;
- spec2=spec2/numrecord;
- ex2time=ex2time/numrecord
- crossspec=crossspec/numrecord;
- transferfn=crossspec./spec2;
- coherence=abs(crossspec).^2./(spec1.*spec2);
- ex2spec=(spec1(1)/2+sum(spec1(2:num/2)))*delfreq %compute the mean square using the power spectrum by integrating over all frequencies.
- %estimate the probability density of the random signal.
- numbin=100;
- bin(1:numbin+1)=0;
- xmax=max(audioarray(:,1));
- xmin=min(audioarray(:,1));
- delx=(xmax-xmin)/numbin;
- xbin=linspace(xmin,xmax,numbin+1);
- for i=1:num*numrecord
- ibin=round((audioarray(i,1)-xmin)/delx)+1;
- bin(ibin)=bin(ibin)+1;
- end
- bin=bin./(delx*num*numrecord);
- sumofprobability=sum(bin)*delx
- figure ('name','time psd and density')
- subplot(3,1,1)
- plot(time,rtime1)
- ylabel('signal')
- xlabel('time in seconds')
- grid on
- subplot(3,1,2)
- semilogy(freq(1:num/2),spec1(1:num/2),freq(1:num/2),spec2(1:num/2))
- ylabel('psd')
- xlabel('frequency in Hz')
- grid on
- subplot(3,1,3)
- plot(xbin,bin)
- ylabel('prob dens')
- xlabel('bin number')
- grid on
- figure('name','transfer function 2/1')
- subplot(2,1,1)
- loglog(freq(1:num/2),abs(transferfn(1:num/2)))
- ylabel('transfer funtion 2/1')
- xlabel('frequency in Hz')
- grid on
- subplot(2,1,2)
- semilogx(freq(1:num/2),angle(transferfn(1:num/2)))
- ylabel('phase of xferfn in rad')
- xlabel('frequency in Hz')
- grid on
- figure('name','coherence')
- semilogx(freq(1:num/2),coherence(1:num/2))
- ylabel('coherence')
- xlabel('frequency (Hz)')
- grid on
- figure('name','Power Spectral Density G(f) (unit^2/Hz)')
- loglog(freq(1:num/2),spec1(1:num/2),freq(1:num/2),spec2(1:num/2))
- ylabel('psd')
- xlabel('frequency in Hz')
- legend('spec1','spec2')
- grid on
- figure('name','psd dB')
- semilogx(freq(1:num/2),10*log10(spec1(1:num/2)),freq(1:num/2),10*log10(spec2(1:num/2)))
- ylabel('psd')
- xlabel('frequency in Hz')
- legend('spec1','spec2')
- grid on
- %Write an ascii file so the data may be plotted in another program.
- %dataout=[freq(1:num/2)',spec1(1:num/2),spec2(1:num/2),coherence(1:num/2),abs(transferfn(1:num/2)),angle(transferfn(1:num/2))*180/pi];
- %save 'soundcard.dat' dataout -ascii
- %%
- % Calculate overall sound pressure level using the program.
- MSP_time = ex2time; MSP_spec = ex2spec; P_ref = 20e-6;
- SPL_time = 10 * log10(MSP_time / P_ref^2)
- SPL_spec = 10 * log10(MSP_spec / P_ref^2)
- % Calculation of 1/3 octave band center frequencies. Start with the lower
- % frequency at band 43
- f(1)=17959.393;
- % Equation found online to calculate all lower frequencies given a starting
- % frequency.
- for n=1:32
- f(n+1)=f(n)/2^(1/3);
- end
- Lower_f(1:32)=f(32:-1:1);
- % Upper frequencies are the same as the lower frequencies except ahead by
- % 1, and the final upper frequency at band 43 is inserted.
- Upper_f(1:32)=[Lower_f(2:32),22627.417];
- % Table of the 1/3 octave band lower and upper frequencies.
- table((12:43)',Lower_f',Upper_f','VariableNames' , {'Band' 'Lower' 'Upper'})
- % Find the index number corresponding to the 'freq' array of the program.
- indexL=ceil(Lower_f ./ fs.*length(freq));
- indexU=ceil(Upper_f ./ fs.*length(freq));
- % Determining the mean square pressure for each 1/3 octave band (for band
- % numbers 12 to 43).
- for i=1:32
- MSP_OB(i)=sum(spec1(indexL(i):indexU(i)))*delfreq;
- end
- % Overall 1/3 octave band sound pressure level.
- OA13_SPL = 10*log10(sum(MSP_OB)/P_ref^2)
- % Percent difference comparison of the two overall sound pressure levels.
- diff = abs((SPL_time - OA13_SPL) / mean([SPL_time , OA13_SPL])) * 100
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement