Narrow To N th Octave

admin

https://ww2.mathworks.cn/matlabcentral/fileexchange/52386-narrowtonthoctave

This function takes narrow band data and converts it to 1/n octave data. Typical uses will be where data were acquired with a constant delta frequency (e.g. from time domain data converted to frequency domain with a FFT), but the user will want this data in octave or 3rd octave form. This function solves the problem of how to convert the data.
[OctaveData,OctaveCenterFrequencies,Flow,Fhigh] = NarrowToNthOctave(narrowFreqArray,narrowdBArray,1)
[thirdOctaveData,ThirdOctaveCenterFrequencies,Flow,Fhigh] = NarrowToNthOctave(narrowFreqArray,narrowdBArray,3)

1. %Program to convert narrowband data to 1/n octaveband data. Just send it an
   
2. %array of band center frequencies, the dB values to sort, and the nth
   
3. %octave band you want. This can be used with any dB values (e.g. Sound
   
4. %pressure level (SPL), sound power (Lw), sound intensity (Li), etc..
   
5. %
   
6. %Important Note 1: The given center frequencies must be a constant df apart
   
7. %
   
8. %Important Note 2: If your data was windowed during the FFT process and you
   
9. %used the amplitude correction factor (ACF), then you need to adjust your
   
10. %data to the energy correction factor when converting to octave spectra. See
    
11. %the first couple lines of the function to adjust this.
    
12. %More reading here: http://blog.prosig.com/2009/09/01/amplitude-and-energy-correction-a-brief-summary/
    
13. %
    
14. %-Examples:
    
15. %[sortedData,Fc,Flow,Fhigh] = NarrowToNthOctave([490,510,500],[60,60,60],1)
    
16. %[sortedData,Fc] = NarrowToNthOctave([490,510,500],[60,60,60],3)
    
17. %[sortedData,Fc] = NarrowToNthOctave(freqArray,dBArray,3)
    
18. %
    
19. %-V4 TMB 12/1/2015
    
20. function [sortedData,Fc,Flow,Fhigh] = NarrowToNthOctave(arrayOfNBCenterFreqs,arrayOfdBToConvert,n)
    
21. df=arrayOfNBCenterFreqs(2)-arrayOfNBCenterFreqs(1);
    
22. %Determine Initial Center frequency (Fc)
    
23. previous=1000*2^(1/n);
    
24. lowFc=1000; %Start at 1000 Hz and find lowest Fc
    
25. while (previous-lowFc) > df
    
26.     previous=lowFc;
    
27.     lowFc=lowFc/2^(1/n);
    
28. end
    
29. %Compute center frequencies
    
30. ii=1;
    
31. Fc(ii)=lowFc;
    
32. while Fc(ii) <= max(arrayOfNBCenterFreqs)
    
33.     ii=ii+1;
    
34.     Fc(ii)=2^(1/n)*Fc(ii-1);
    
35. end
    
36. %Compute high, low frequencies, and bandwidth (BW) from center frequencies
    
37.     Flow=Fc/sqrt(2^(1/n));
    
38.     Fhigh=Fc*sqrt(2^(1/n));
    
39.     BW=Fhigh-Flow;
    
40. %Sort data in frequency bins
    
41. sortedData=zeros(1,length(Fc)); zeroFreqValue=0;
    
42. for jj=1:length(arrayOfdBToConvert) %Do for each value given
    
43.     for kk=1:length(Fc) %Check each Fc bin to see where value should be placed
    
44.         if arrayOfNBCenterFreqs(jj)>=Flow(kk) && arrayOfNBCenterFreqs(jj)<Fhigh(kk) %Find place. In the rare case a given center freq = flow(ii)/fhigh(ii+1) then sum the value in flow
    
45.             if sortedData(kk)==0 %if no values has been added to the band then set initial value
    
46.                 sortedData(kk)=arrayOfdBToConvert(jj);
    
47.                 N(kk)=1;
    
48.             else %else sum in the value
    
49.                 sortedData(kk)=10*log10(sum(10.^([sortedData(kk) arrayOfdBToConvert(jj)]./10))); %dB add value to that bin
    
50.                 N(kk)=N(kk)+1;
    
51.             end
    
52.         else %else check to see if the value belongs in the 0 Hz band
    
53.             if (kk==1) && arrayOfNBCenterFreqs(jj)<Flow(kk)
    
54.                 if zeroFreqValue==0
    
55.                     zeroFreqValue=arrayOfdBToConvert(jj);
    
56.                 else
    
57.                     zeroFreqValue=10*log10(sum(10.^([zeroFreqValue arrayOfdBToConvert(jj)]./10))); %dB add value to that bin
    
58.                 end
    
59.             end
    
60.         end
    
61.     end
    
62. end
    
63. %Apply bin correction factor for when center frequencies are not integer
    
64. %multiples of the bandwidth (BW)
    
65. for ll=1:length(N) %Look through each bin that had values added to it
    
66.     if (N(ll)~=0) && (floor(BW(ll)/df)<=N(ll)) %If values were added and there are enough values added to make applying the correction appropriate then apply correction
    
67.         sortedData(ll)=sortedData(ll)+10*log10(BW(ll)/(df*N(ll)));
    
68.     end
    
69. end
    
70. %If DC value was given, add to sorted array
    
71. if zeroFreqValue~=0
    
72.     sortedData=[zeroFreqValue sortedData];
    
73.     Fc=[0 Fc];
    
74.     Flow=[0 Flow];
    
75.     Fhigh=[0 Fhigh];
    
76. end
    
77. %Remove unused bins at beginning if only higher freqs were given
    
78. removeThisMany=0;
    
79. for ll=1:length(sortedData)
    
80.     if sortedData(ll)==0;
    
81.         removeThisMany=1+removeThisMany;
    
82.     else
    
83.         break
    
84.     end
    
85. end
    
86. sortedData=sortedData(removeThisMany+1:end);
    
87. Fc=Fc(removeThisMany+1:end);
    
88. Flow=Flow(removeThisMany+1:end);
    
89. Fhigh=Fhigh(removeThisMany+1:end);
    
90. %Remove bin at end if unused
    
91. if sortedData(end)==0
    
92.     sortedData=sortedData(1:end-1);
    
93.     Fc=Fc(1:end-1);
    
94.     Flow=Flow(1:end-1);
    
95.     Fhigh=Fhigh(1:end-1);
    
96. end
    
97. end

复制代码