%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MatLab script "ANSS_noise_rms_rev4.m":
%
% Sample MatLab code for ANSS-recommended sensor self-noise
% rms calculation in the case where ambient noise is known
% to be well below self noise.
%
% Use for low-gain accelerometers and ambient ground motion
% << sensor noise. (Use 2-channel or 3-channel method when
% ambient ground motion >~ sensor noise.)
%
% Input file defined in InputFname [any single column of
% acceleration values, suggested minimum 30 minutes and
% 160000 points of data at 200 sps in order to resolve
% low frequencies and optimize ensemble averaging].
%
% CAVEAT: This script rounds down the input file duration
% to the nearest modulo(nSegs) seconds so that the nSegs
% windows will consist of a round number of seconds in
% duration. nSegs = 10 is typical, yielding rounding to
% the nearest 10-s mark below the original trace length.
%
% Output is three plots plus file 'accel_output.txt' [as four
% text columns, comma-delimited: Frequency, PSD (cm/s/s/Hz),
% 1/2-octave rms (cm/s/s), 1/3-octave rms (cm/s/s) ]. Only the
% rms and PSD are ANSS recognized standards.
%
% Calculations use Welch's Method to estimate an averaged PSD,
% then multiply by bandwidth (1/3 or 1/2-octave) and take
% square root to get rms estimate.
%
% USAGE:
% Either cut and paste the contents of this file to MatLab,
% or in MatLab, select the directory containing this script,
% then issue "ANSS_noise_rms_rev4" as a command (no quotes).
% For example:
%
% close all;
% clear all;
% cd ('/Users/jrevans/NSMP/TIC_WGD8/Analysis/RMS Computation');
% ANSS_noise_rms_rev4
%
% Calls: (only standard MatLab routines)
%
% History:
% Rev. 0: R. Nigbor, UCLA, 24 June 2006
% Rev. 1: JREvans, USGS, Menlo Park, 26 June 2006
% Rev. 2: R. Nigbor, UCLA, 26 June 2006
% Rev. 3: JREvans, USGS, Menlo Park, 28 June 2006
% Rev. 4: JREvans, USGS, Menlo Park, 13 Nov 2006
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diary off
close all;
clear all;
fclose('all');
CR = sprintf('\n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% User parameters (adjust as needed/desired)
ActiveDir = '/Users/jrevans/TREMOR_SOS/NetQuakes/OFR_SaPGVetc/SelfNoise';
cd (ActiveDir); % Directory for input and output
InputFname = 'accel_input.txt'; % Input-file name
OutputFname = 'accel_output.csv';% Output comma-separated-file name
% Plots will be output as PostScript
% with names based on InputFname.
one_g = 980.665; % (NB: Ideally, this value should
% be set to the local gravitational
% acceleration at the site where the
% accelerometer was calibrated, unless
% that calibration has been corrected
% to this international standard "g".
%
% However, the differences are only
% parts per thousand and not likely
% to be significant.)
SampleRate = 200; % Sample rate in samples per second.
MaxFreq = SampleRate / 2.; % Maximum frequency for plots, Nyquist
% unless changed.
ScaleToGal = 1.0; % Multiplicative factor to scale data
% to cm/s/s (set equal to one_g if
% input data are in g's)
doSmooth = 1; % For any nonzero value, smoothes the
% rms partial-octave spectra by running
% means of five. Full RMS spectrum is
% NOT smoothed.
SavePlot = 1; % For any nonzero value, saves the plots
% to Postscript files in ActiveDir.
ClipLevel = 3.5; % Sensor clip level in g's.
ClipLevel = ClipLevel * one_g; % Sensor clip level in cm/s/s.
nSegs = 10; % Number of windows to subset the trace
% into (disregarding overlap). We
% suggest at least 10 and round numbers.
% Plot-title lines 1 through 3; change the SECOND one appropriately.
pagetitle1 = 'ANSS-Recommended Low-Gain Accelerometer Noise Analysis';
pagetitle2 = [ 'Sensor Model MMMMM, Serial Number NNNNN; ' ...
'Test Date dd/mmm/yyyy' ];
if (doSmooth)
pagetitle3 = [ '(Smoothed; ' ...
strrep([ 'File "' ActiveDir '/' InputFname '")' ], ...
'_', ' ') ];
else
pagetitle3 = [ '(UNsmoothed; ' ...
strrep([ 'File "' ActiveDir '/' InputFname '")' ], ...
'_', ' ') ];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Read the data set and trim off the tail end points to a size
% such that subsegments are an integer number of seconds long.
[ 'Reading acceleration trace from file' CR ...
' ' ActiveDir '/' InputFname CR ...
' (this may take a few minutes) ...' ]
Accel = load(InputFname); % A column vector
[ nPtsOld, nCol ] = size(Accel);
if (nCol ~= 1)
error('Acceleration data must be column vector.');
end
% Trim to a duration that will result in pwelch() windows
% with whole-second durations.
nSec = (nPtsOld / SampleRate) - mod(nPtsOld / SampleRate, nSegs);
nPts = nSec * SampleRate;
Accel = Accel(1:nPts);
% Notify of trimming
if (nPts < nPtsOld)
disp([ 'To fit subsegments evenly, trimmed off end of time series by ' ...
num2str(nPtsOld - nPts) ...
' points' CR ...
' yielding ' num2str(nPts) ' points (' ...
num2str(nPts / SampleRate) ' seconds; ' ...
num2str(nPts / SampleRate / 60) ' minutes).' ]);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Additional constants calculated automatically or set directly
nFFT = nPts / nSegs; % FFT block size for ensemble
% nFFT = 5000; % averaging in pwelch()
nOverlap = nFFT / 2; % Overlap size for ensemble
% nOverlap = 2500; % averaging in pwelch()
Nyquist = SampleRate / 2;
LowBin = (Nyquist) / (nFFT / 2);
% Summarize computational arrangements
if (rem(nOverlap, 1) ~= 0)
error([ CR 'Sample rates that do not yield even-length ' ...
'windows' CR ' when the time series is rounded ' ...
'down to the nearest ' CR ' ' num2str(nSegs) ...
' seconds are not allowed.' CR ]);
end
disp([ 'In pwelch(), using windows of ' num2str(nFFT) ' points (' ...
num2str(nFFT / SampleRate) ' s)' CR ' overlapping by ' ...
num2str(nOverlap) ' points (' ...
num2str(nOverlap / SampleRate) ' s).' ]);
disp([ ' Lowest frequency bin = ' num2str(LowBin) ...
' Hz; Nyquist = ' num2str(Nyquist) ' Hz.' ]);
if ((LowBin) > 0.01)
disp([ 'Caution: lowest AC frequency bin is above 0.01 Hz.' ]);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Correct Accel as needed
% Remove mean, and scale to cm/s/s
Accel = detrend(Accel, 'constant'); % remove mean
Accel = Accel .* ScaleToGal; % scale to cm/s/s
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate PSD using Welch's Method ('help pwelch()' for info)
% Use Hanning window the same length as each ensemble
disp([ CR 'For adequately long noise samples, this may take a minute or so ...' ]);
[ PSD_Acc, Freq ] = ...
pwelch(Accel, hanning(nFFT, 'periodic'), nOverlap, nFFT, ...
SampleRate);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate time and frequency domain rms for comparison check
rms_t = sqrt(mean(Accel .^ 2));
rms_f = sqrt(trapz(PSD_Acc) * SampleRate / nFFT);
disp([ CR 'Full-band rms in the frequency domain is ' num2str(rms_f) ...
' cm/s/s;' CR ' rms in the time domain is ' ...
num2str(rms_t) ' cm/s/s;' CR ' so the percentage ' ...
'error, 100*(frequency/time - 1), is ' ...
num2str(100 * ((rms_f/rms_t) - 1), 2) '%' CR ]);
% Calculate rms arrays
ThirdOctave = sqrt(PSD_Acc .* Freq * (2^(1/3) - 1)); % 1/3-octave rms
HalfOctave = sqrt(PSD_Acc .* Freq * (2^(1/2) - 1)); % 1/2-octave rms
RMS = sqrt(PSD_Acc); % rms spectrum
% Smooth partial-octave spectra using a 5-point running mean (optional)
if (doSmooth)
ThirdOctave = smooth(ThirdOctave);
HalfOctave = smooth(HalfOctave);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Plots
% Make time vector for plotting
Tmax = (nPts - 1) / SampleRate;
Time = linspace(0, Tmax, nPts);
fmax = 10 ^ floor(log10(MaxFreq));
fmin = 10 ^ floor(log10(Freq(2)));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1);
scrsz = get(0,'ScreenSize');
adjx = round(0.125 * scrsz(3));
adjy = round(0.100 * scrsz(4));
scrsz = [ adjx adjy (scrsz(3) - (2 * adjx)) (scrsz(4) - (2 * adjy)) ];
% Fill screen with plot (may not work well on PC -- comment it out)
set(gcf, 'Position', scrsz);
subplot(2, 1, 1);
amax = max(abs(Accel));
plot(Time, Accel);
axis([ 0 Tmax -amax amax ]);
grid on;
title({pagetitle1;pagetitle2;pagetitle3}, 'FontSize', 14, ...
'FontWeight', 'bold');
ylabel('Acceleration cm/s^2', 'FontSize', 12, 'FontWeight', 'bold');
xlabel('Time (s)', 'FontSize', 12, 'FontWeight', 'bold');
legend(strcat('Time Domain rms = ', num2str(rms_t, 5)), ...
'Location', 'SouthWest')
subplot(2, 1, 2);
NonzeroPSD = nonzeros(PSD_Acc);
pmax = 10 ^ ceil(log10(max(NonzeroPSD)));
pmin = 10 ^ floor(log10(min(NonzeroPSD)));
loglog(Freq, PSD_Acc, 'LineWidth', 1.25);
axis([ fmin fmax pmin pmax ]);
grid on
title('Power Spectral Density (Welch''s Method with Hanning Windows)', ...
'FontSize', 14, 'FontWeight', 'bold');
ylabel('PSD ((cm/s^2)^2/Hz)', 'FontSize', 12, 'FontWeight', 'bold');
xlabel('Frequency (Hz)', 'FontSize', 12, 'FontWeight', 'bold');
legend(strcat('Frequency Domain rms = ', num2str(rms_f, 5)), ...
'Location', 'SouthWest')
if (SavePlot)
orient landscape, print('-dpsc', ...
[ strrep(InputFname, '.txt', '') '_Fig_1.ps' ]);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2);
scrsz = get(0,'ScreenSize');
adjx = round(0.125 * scrsz(3));
adjy = round(0.100 * scrsz(4));
scrsz = [ adjx adjy (scrsz(3) - (2 * adjx)) (scrsz(4) - (2 * adjy)) ];
% Fill screen with plot (may not work well on PC -- comment it out)
set(gcf, 'Position', scrsz);
Fclip = [ fmin fmax ];
Aclip = [ (ClipLevel * sqrt(0.5)) (ClipLevel * sqrt(0.5)) ];
NonzeroAcc = nonzeros([ ThirdOctave' HalfOctave' ...
(ClipLevel * sqrt(0.5)) ])';
omax = 10 ^ ceil(log10(max(NonzeroAcc)));
omin = 10 ^ floor(log10(min(NonzeroAcc)));
loglog(Freq, ThirdOctave, 'c', Freq, HalfOctave, 'b',...
Fclip, Aclip, 'm', 'LineWidth', 1.25);
axis([ fmin fmax omin omax ]);
grid on;
title({pagetitle1;pagetitle2;pagetitle3}, ...
'FontSize', 14, 'FontWeight', 'bold');
ylabel('Narrowband rms Acceleration (cm/s^2)', ...
'FontSize', 12, 'FontWeight', 'bold');
xlabel('Frequency (Hz)', 'FontSize', 12, ...
'FontWeight', 'bold');
legend('1/3-Octave rms', '1/2-Octave rms', ...
'Clip Level, rms Sine', 'Location', 'West');
if (SavePlot)
orient landscape, print('-dpsc', ...
[strrep(InputFname, '.txt', '') '_Fig_2.ps']);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(3);
scrsz = get(0,'ScreenSize');
adjx = round(0.125 * scrsz(3));
adjy = round(0.100 * scrsz(4));
scrsz = [ adjx adjy (scrsz(3) - (2 * adjx)) (scrsz(4) - (2 * adjy)) ];
% Fill screen with plot (may not work well on PC -- comment it out)
set(gcf, 'Position', scrsz);
NonzeroRMS = nonzeros([ RMS' (ClipLevel * sqrt(0.5)) ])';
omax = 10 ^ ceil(log10(max(NonzeroRMS)));
omin = 10 ^ floor(log10(min(NonzeroRMS)));
loglog(Freq, RMS, 'r', Fclip, Aclip, 'm', 'LineWidth', 1.25);
axis([ fmin fmax omin omax ]);
grid on;
% Raw RMS spectrum is NEVER smoothed
pagetitle3x = [ '(UNsmoothed; ' strrep([ 'File "' ActiveDir '/' InputFname '")' ], '_', ' ') ];
title({pagetitle1;pagetitle2;pagetitle3x}, 'FontSize', 14, ...
'FontWeight', 'bold');
ylabel('Per-bin rms Acceleration (cm/s^2)/sqrt(Hz)', ...
'FontSize', 12, 'FontWeight', 'bold');
xlabel('Frequency (Hz)', 'FontSize', 12, 'FontWeight', 'bold');
legend('rms', 'Clip Level, rms Sine', 'Location', 'SouthWest');
if (SavePlot)
orient landscape, print('-dpsc', ...
[strrep(InputFname, '.txt', '') '_Fig_3.ps']);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Output to comma-delimited file, Freq, PSD, 1/3-octave, 1/2-octave
dlmwrite(OutputFname, [ Freq, PSD_Acc, RMS, ThirdOctave, HalfOctave ]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%