Calibrated Spectral Analysis: [SPa]=test_pressure_spectra(demos, Fs, fc, bin_size, num_averages, win_type, flag1 ) - File Exchange

Code covered by the BSD License

Highlights from
Calibrated Spectral Analysis

ACdsgn(Fs)
ArgStruct=parseArgs(args,ArgS... Helper function for parsing varargin.
[Fsn, p, q, errors]=get_p_q2(...
[SP, f, bin_size, num_average... % pressure_spectra: Calculates an accurate estimate of the pressure spectra
[SP, f, num_averages_out]=spe... % spectra_estimate: Is a rough estimate of the pressure spectra
[SP2, mean_array2]=sub_mean(S... % sub_mean: Removes the running average from a time record given a sampling rate and high pass cutoff frequency.
[SPa]=test_pressure_spectra(d... % test_spectra_estimate: runs demos for the pressure spectra program.
[bin_size, num_averages_out, ... % number_of_averages: Calculates the number of points not overlapped from the array size, bin size, and number of averages
[bz, az]=bessel_digital(Fs, F... % bessel_digital: creates a digital low pass bessel filter of order n
[cfa, SP, f, f_cal]=mic_calib... % mic_calib: Uses a flat top window to calibrate using A-weighted or Linear weighting
[ftwcf]=window_correction_fac... % window_correction_factor: Computes the factor for calibrating a Fourier Transform given specific processing parameters
[m2]=geospace(a, b, n, flag) % geospace: caculates a geometric sequence or psuedogeometric sequence from a to b with n elements
[prms]=rms_val(p, dim) % rms_val: Calculates the rms value along a specific dimension
[varargout]=convert_double(va... % This program converts the inputs into double precision arrays. Then
[w]=flat_top(N, type) % Flat top windows are used for calibration, because the wide main lobe
[y, x, a]=match_height_and_sl... % match_height_and_slopes2: creates a quartic with specifed height and slope at the end points.
[y2, num_settle_pts, settling... % filter_settling_data: Creates data to append to a time record for settling a filter
[y2]=remove_filter_settling_d... % remove_filter_settling_data: removes data added to time records to settle the filter
[yAC, errors]=ACweight_time_f... % ACweight_time_filter: Applies an A or C weighting time filter to a sound recording
[y]=moving(x,m,fun) MOVING will compute moving averages of order n (best taken as odd)
[y_out, b, a]=bessel_antialia... % bessel_antialias: applies an antialiasing digital Bessel filter
[y_out, t_out, b, a]=bessel_d... % bessel_down_sample: applies an antialiasing digital Bessel filter
[y_out, x_out, y_in]=resample... % resample_interp3: resamples using interp1 with additional options
h=subaxis(varargin) SUBAXIS Create axes in tiled positions. (just like subplot)
loc=LMSloc(X)
View all files

from Calibrated Spectral Analysis by Edward Zechmann
Simple Fourier Spectral Analysis of sound pressure time record.

[SPa]=test_pressure_spectra(demos, Fs, fc, bin_size, num_averages, win_type, flag1 )

function [SPa]=test_pressure_spectra(demos, Fs, fc, bin_size, num_averages, win_type, flag1 )
% % test_spectra_estimate: runs demos for the pressure spectra program.
% %
% % Syntax:
% %
% % [SPa]=test_pressure_spectra(demos, Fs, fc, bin_size, num_averages, win_type, flag1);
% %
% % ***********************************************************
% %
% % Description
% %
% % test_pressure_spectra creates a siusoid then calculates and plots the
% % pressure spectra.  This allows the user to see teh pressur espectra for
% % different sampling rates center frequencies, bin_sizes, 
% % number of averages, and window types.  
% % 
% % Three different demos can be run by setting demos to 1, 2, or 3.    
% % Assumes a test sinusoid signal level of 114 dB.
% % 
% % 
% % 
% % 
% % ***********************************************************
% %
% % Input Variables
% %
% % demos=1;    % 250 Hz 114 dB test sinusoid with 
% %             % Fs=44100;
% %             % N=65536;
% %             % bin_size=65536;
% %             % num_averages=1;
% %             % blackman window
% %
% % demos=2;    % 250 Hz 114 dB test sinusoid with
% %             % Fs=44100;
% %             % N=16*65536;
% %             % bin_size=65536;
% %             % num_averages=100;
% %             % hanning window
% %
% % demos=3;    % 250 Hz 114 dB test sinusoid with
% %             % Fs=44100;
% %             % N=16*65536;
% %             % bin_size=65536;
% %             % num_averages=10000;
% %             % hanning window
% %
% %
% % Fs (Hz) is the Sampling Rate.
% %
% % fc (Hz) is the center frequency for a sinusoidal test signal.
% %
% % bin_size is the number_of points in each fft should be divisible by 2.
% %
% % num_averages is the Number of time averages.
% %
% % win_type is the type of window for tapering the beginning and ending
% %             of the time records to zero before computings the FFTs.
% %
% %      List of supported windows other windows may work too.
% %
% %           blackman
% %           chebwin
% %           flat_top
% %           flattopwin
% %           gausswin
% %           hamming
% %           hann
% %           hanning
% %           kaiser
% %           tukeywin
% %
% %
% % flag1=1 forces the progrmam to calculate the maximum number of averages
% %             using an overlap of one data point.
% %
% %
% % ***********************************************************
% %
% % Output Variables
% %
% % SPa (Pa) is the rms sound pressure.
% %
% %
% %
% % ***********************************************************
% %
% %
% %
% % List of Dependent Subprograms for
% % test_pressure_spectra
% %
% % FEX ID# is the File ID on the Matlab Central File Exchange
% %
% %
% % Program Name   Author   FEX ID#
% %  1) ACdsgn		Edward L. Zechmann
% %  2) ACweight_time_filter		Edward L. Zechmann
% %  3) bessel_antialias		Edward L. Zechmann
% %  4) bessel_digital		Edward L. Zechmann
% %  5) bessel_down_sample		Edward L. Zechmann
% %  6) convert_double		Edward L. Zechmann
% %  7) filter_settling_data3		Edward L. Zechmann
% %  8) flat_top		Edward L. Zechmann
% %  9) geospace		Edward L. Zechmann
% % 10) get_p_q2		Edward L. Zechmann
% % 11) LMSloc		Alexandros Leontitsis		801
% % 12) match_height_and_slopes2		Edward L. Zechmann
% % 13) moving		Aslak Grinsted		8251
% % 14) number_of_averages		Edward L. Zechmann
% % 15) pressure_spectra		Edward L. Zechmann
% % 16) remove_filter_settling_data		Edward L. Zechmann
% % 17) resample_interp3		Edward L. Zechmann
% % 18) rms_val		Edward L. Zechmann
% % 19) spectra_estimate		Edward L. Zechmann
% % 20) sub_mean		Edward L. Zechmann
% % 21) window_correction_factor		Edward L. Zechmann
% %
% %
% %
% % ***********************************************************
% %
% % This program was written by Edward L. Zechmann
% %
% %     date  28 February   2008
% %
% % modified  29 September  2008  Added a dependent function list
% %
% % modified  19 February   2011    Removed subtracting off the running
% %                                 average
% %
% %
% % ***********************************************************
% %
% % Feel free to modify this code.
% %
% % See Also: spectra_estimate, fft, Aweight_time_filter,
% %           Cweight_time_filter, AC_weight_NB
% %
% %
% %
% %

wintype1={'blackman', 'chebwin', 'flat_top', 'flattopwin', 'gausswin', 'hamming', 'hann', 'kaiser', 'tukeywin'};
m1=length(wintype1);

if nargin >= 1 && ~isempty(demos) && isnumeric(demos) && (demos >= 1)  && (demos <= 3)
    
    demos=round(demos);
    
    switch demos
        
        case 1
            
            Fs=44100;
            N=65536;
            bin_size=65536;
            num_averages=1;
            f_cal=250;
            flag1=0;
            amp=0.00002*10^(114/20)*sqrt(2);
            x=amp*sin(f_cal*2*pi/Fs*(1:N));
            m1=length(wintype1);
            
        case 2
            
            Fs=44100;
            N=16*65536;
            bin_size=65536;
            num_averages=100;
            f_cal=250;
            flag1=0;
            amp=0.00002*10^(114/20)*sqrt(2);
            m1=1;
            wintype1={'hann'};
            x=amp*randn(1, N);
            
            %x=[x; randn(1, N); y];
        case 3
            
            Fs=44100;
            N=16*65536;
            bin_size=65536;
            num_averages=10000;
            f_cal=250;
            flag1=0;
            amp=0.00002*10^(114/20)*sqrt(2);
            m1=1;
            wintype1={'hann'};
            x=amp*randn(1, N);
            
            
        otherwise
            
            Fs=44100;
            N=65536;
            bin_size=65536;
            num_averages=1;
            f_cal=250;
            flag1=0;
            amp=0.00002*10^(114/20)*sqrt(2);
            x=amp*sin(f_cal*2*pi/Fs*(1:N));
            m1=length(wintype1);
    end
    
else
    
    if nargin < 1 || isempty(demos) || ~isnumeric(demos)
        demos=[];
    end
    
    if nargin < 2 || isempty(Fs) || ~isnumeric(Fs)
        Fs=50000;
    end
    
    if nargin < 3 || isempty(Fc) || ~isnumeric(Fc)
        fc=250;
    end
    
    if nargin < 4 || isempty(bin_size) || ~isnumeric(bin_size)
        bin_size=length(x);
    end
    
    if nargin < 5 || isempty(num_averages) || ~isnumeric(num_averages)
        num_averages=1;
    end
    
    if nargin < 6 || isempty(win_type) || ~ischar(win_type)
        win_type='hann';
    end
    
    if nargin < 7 || isempty(flag1) || ~isnumeric(flag1)
        flag1=0;
    end
    
    
    
    m1=1; % number of windows to apply
    amp=0.00002*10^(114/20)*sqrt(2);
    x=amp*sin(fc*2*pi/Fs*(1:N));
    
end






[num_channels, num_data_points]=size(x);


for e1=1:m1;
    
    [SP, f]=pressure_spectra(x, Fs, bin_size, num_averages, wintype1{e1}, flag1 );
    
    if isequal(e1,1)
        SPa=zeros( m1, num_channels, length(SP));
    end
    
    for e2=1:num_channels;
        for e3=1:length(SP);
            SPa(e1, e2, e3)=20*log10(1/0.00002*SP(e2, e3));
        end
    end
    
end


% plot with subaxis
sh=0.0;
sv=0.0;
ml=0.14;
mr=0.1;
mt=0.08;
mb=0.12;



l_min=0.1*min(min(floor(10*(SPa-10))));
l_max=0.1*max(max(ceil(10*(SPa+10))));

for e1 =1:m1;
    
    for e2=1:num_channels;
        for e3=1:length(SP);
            buf(e2, e3)=SPa(e1, e2, e3);
        end
    end
    
    % plot with subaxis
    subaxis(m1, 1, e1, 'sh', sh, 'sv', sv , 'pl', 0, 'pr', 0, 'pt', 0, 'pb', 0, 'ml', ml, 'mr', mr, 'mt', mt, 'mb', mb);
    
    semilogx(f, buf, 'LineWidth', 2);
    
    if e1==1
        ylabel('Sound (dB ref. 20\muPa)', 'Fontsize', 16);
    end
    
    xlabel('Frequency (Hz)', 'Fontsize', 16);
    title( ['Test Spectra for ', wintype1{e1}, ' Window'], 'Interpreter', 'none', 'Fontsize', 20 );
    set(gca, 'Fontsize', 16);
    text(0.8*max(f), l_max-0.02*(l_max-l_min), ['Max Level ', num2str(0.01*max(max(round(100*SPa(e1, :, :)))))], 'color', [0 0 0], 'Fontsize', 16, 'HorizontalAlignment', 'right', 'VerticalAlignment', 'top');
    
    ylim([l_min l_max]);
    xlim([min(f) max(f)]);
    
end

Highlights from Calibrated Spectral Analysis

Highlights from
Calibrated Spectral Analysis