Advertisement
Guest User

Untitled

a guest
Sep 17th, 2014
213
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 6.54 KB | None | 0 0
  1. % Modified by: Brian Cai
  2. % ME 622: Advanced Acoustics
  3. % Assignment 2: Problem 2
  4. % Due September 18, 2014
  5.  
  6. % This m file computes the autospectrum, prob density, and transfer
  7. % funtions from a time domain signal acquired by the computer. It is
  8. % written using Matlab Version 7.0. Note that to acquire data using this
  9. % progam the user must make sure that the desired recording device is
  10. % "selected" from the "Record Control" dialog box in Windows. To do this,
  11. % right click on the loudspeaker icon in the quick launch section at the
  12. % right end of the task bar and select "Open Volume Control." Click
  13. % "options" and select "properties" and then select "recording". Then
  14. % select "Line-In" or "Microphone" or "CD Audio" etc. depending on which
  15. % input you want to record from. If you don't have the loudspeaker icon in
  16. % your quick launch section you can select the recording devices from
  17. % Control Panel>Sounds and Audio Devices.
  18. % R. N. Miles
  19. clear all; close all; clc
  20.  
  21. % The number of records to acquire.
  22. numrecord = 50;
  23.  
  24. % num is the total number of time points per record. This should be a power
  25. % of 2 for the FFT algorithm.
  26. num = 16384;
  27.  
  28. % The sample rate in samples/second. This is a reasonable rate for audio
  29. % signals.
  30. fs = 44100;
  31. % fs = 96000;
  32.  
  33. % The amount of time in each record.
  34. Trecord = num / fs;
  35.  
  36. % Total time in seconds.
  37. totaltime = numrecord * Trecord;
  38.  
  39. % Acquire data from the sound card at the rate fs (samples/s) using 16
  40. % bits, two channels.
  41. recorder = audiorecorder(fs , 16 , 2);
  42. % recorder = audiorecorder(fs , 24 , 2);
  43. recordblocking(recorder , totaltime);
  44. audioarray = getaudiodata(recorder);
  45.  
  46. % The rest of the code remains unchanged.
  47. delfreq=1/Trecord; %The frequency spacing in Hertz.
  48. numt=length(audioarray(:,1));
  49. time=linspace(0,Trecord,num); %Time for first record
  50. figure ('name','time domain plot')
  51. plot(time,audioarray(1:num,1),time,audioarray(1:num,2))
  52. xlabel('Time in seconds')
  53. ylabel('original signal')
  54. title('Record 1 of original signal')
  55. grid on
  56. delt=1/fs;
  57. freq=linspace(0,fs,num);
  58. spec1=0*freq'; %Create an array of zeros that is the same size as freq.
  59. spec2=spec1;
  60. crossspec=spec1; %Create an array of zeros that is the same size as freq.
  61. ex2time=0;
  62. for irec=1:numrecord
  63.  
  64. %rtime1=audioarray(1+num*(numrecord-1):num+num*(numrecord-1),1); %Store the first num points from channel 1 in an array
  65. %rtime2=audioarray(1+num*(numrecord-1):num+num*(numrecord-1),2); %Store the first num points from channel 2 in an array
  66. rtime1=audioarray(1+num*(irec-1):num+num*(irec-1),1); %Store the first num points from channel 1 in an array
  67. rtime2=audioarray(1+num*(irec-1):num+num*(irec-1),2); %Store the first num points from channel 2 in an array
  68. fft1=fft(rtime1,num);
  69. fft2=fft(rtime2,num);
  70. ex2time=ex2time+sum(rtime1.^2)/num; %Mean square of each time record
  71. spec1=spec1+abs(fft1).^2*2*delt/num; %calculate the autospectrum of the random signal.
  72. spec2=spec2+abs(fft2).^2*2*delt/num; %calculate the autospectrum of the random signal.
  73. % This will be a single-sided spectrum in units^2/Hz
  74. crossspec=crossspec+fft1.*conj(fft2)*2*delt/num;
  75. end
  76. spec1=spec1/numrecord;
  77. spec2=spec2/numrecord;
  78. ex2time=ex2time/numrecord
  79. crossspec=crossspec/numrecord;
  80. transferfn=crossspec./spec2;
  81. coherence=abs(crossspec).^2./(spec1.*spec2);
  82. ex2spec=(spec1(1)/2+sum(spec1(2:num/2)))*delfreq %compute the mean square using the power spectrum by integrating over all frequencies.
  83. %estimate the probability density of the random signal.
  84. numbin=100;
  85. bin(1:numbin+1)=0;
  86. xmax=max(audioarray(:,1));
  87. xmin=min(audioarray(:,1));
  88. delx=(xmax-xmin)/numbin;
  89. xbin=linspace(xmin,xmax,numbin+1);
  90. for i=1:num*numrecord
  91. ibin=round((audioarray(i,1)-xmin)/delx)+1;
  92. bin(ibin)=bin(ibin)+1;
  93. end
  94. bin=bin./(delx*num*numrecord);
  95. sumofprobability=sum(bin)*delx
  96. figure ('name','time psd and density')
  97. subplot(3,1,1)
  98. plot(time,rtime1)
  99. ylabel('signal')
  100. xlabel('time in seconds')
  101. grid on
  102. subplot(3,1,2)
  103. semilogy(freq(1:num/2),spec1(1:num/2),freq(1:num/2),spec2(1:num/2))
  104. ylabel('psd')
  105. xlabel('frequency in Hz')
  106. grid on
  107. subplot(3,1,3)
  108. plot(xbin,bin)
  109. ylabel('prob dens')
  110. xlabel('bin number')
  111. grid on
  112. figure('name','transfer function 2/1')
  113. subplot(2,1,1)
  114. loglog(freq(1:num/2),abs(transferfn(1:num/2)))
  115. ylabel('transfer funtion 2/1')
  116. xlabel('frequency in Hz')
  117. grid on
  118. subplot(2,1,2)
  119. semilogx(freq(1:num/2),angle(transferfn(1:num/2)))
  120. ylabel('phase of xferfn in rad')
  121. xlabel('frequency in Hz')
  122. grid on
  123. figure('name','coherence')
  124. semilogx(freq(1:num/2),coherence(1:num/2))
  125. ylabel('coherence')
  126. xlabel('frequency (Hz)')
  127. grid on
  128. figure('name','Power Spectral Density G(f) (unit^2/Hz)')
  129. loglog(freq(1:num/2),spec1(1:num/2),freq(1:num/2),spec2(1:num/2))
  130. ylabel('psd')
  131. xlabel('frequency in Hz')
  132. legend('spec1','spec2')
  133. grid on
  134. figure('name','psd dB')
  135. semilogx(freq(1:num/2),10*log10(spec1(1:num/2)),freq(1:num/2),10*log10(spec2(1:num/2)))
  136. ylabel('psd')
  137. xlabel('frequency in Hz')
  138. legend('spec1','spec2')
  139. grid on
  140. %Write an ascii file so the data may be plotted in another program.
  141. %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];
  142. %save 'soundcard.dat' dataout -ascii
  143.  
  144. %%
  145.  
  146. % Calculate overall sound pressure level using the program.
  147. MSP_time = ex2time; MSP_spec = ex2spec; P_ref = 20e-6;
  148. SPL_time = 10 * log10(MSP_time / P_ref^2)
  149. SPL_spec = 10 * log10(MSP_spec / P_ref^2)
  150.  
  151. % Calculation of 1/3 octave band center frequencies. Start with the lower
  152. % frequency at band 43
  153. f(1)=17959.393;
  154.  
  155. % Equation found online to calculate all lower frequencies given a starting
  156. % frequency.
  157. for n=1:32
  158. f(n+1)=f(n)/2^(1/3);
  159. end
  160. Lower_f(1:32)=f(32:-1:1);
  161.  
  162. % Upper frequencies are the same as the lower frequencies except ahead by
  163. % 1, and the final upper frequency at band 43 is inserted.
  164. Upper_f(1:32)=[Lower_f(2:32),22627.417];
  165.  
  166. % Table of the 1/3 octave band lower and upper frequencies.
  167. table((12:43)',Lower_f',Upper_f','VariableNames' , {'Band' 'Lower' 'Upper'})
  168.  
  169. % Find the index number corresponding to the 'freq' array of the program.
  170. indexL=ceil(Lower_f ./ fs.*length(freq));
  171. indexU=ceil(Upper_f ./ fs.*length(freq));
  172.  
  173. % Determining the mean square pressure for each 1/3 octave band (for band
  174. % numbers 12 to 43).
  175. for i=1:32
  176. MSP_OB(i)=sum(spec1(indexL(i):indexU(i)))*delfreq;
  177. end
  178.  
  179. % Overall 1/3 octave band sound pressure level.
  180. OA13_SPL = 10*log10(sum(MSP_OB)/P_ref^2)
  181.  
  182. % Percent difference comparison of the two overall sound pressure levels.
  183. diff = abs((SPL_time - OA13_SPL) / mean([SPL_time , OA13_SPL])) * 100
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement