clear all
close all

%Input actual files size used on deployment
samplesPerFile = 20.0 * 60 * 48000.0;   %20 minutes/file * 60 sec/min * 48000 samp/sec

Fs = 48000; %samples per second
dT = 1/Fs;
Vref = 1.5;
LSB = Vref / (2^15 - 1);

dirName = 'N:\DeploymentOct12-09\R2009-285-08-00-26\R00000\';
startFileNum = 80;
numFilesToProcess = 2;
endFileNum = startFileNum + numFilesToProcess - 1;

tStart = tic;

for fileNum=startFileNum:endFileNum
    fileName = strcat(dirName, 'D',  sprintf('%06d', fileNum));
    fid = fopen(fileName);
    fprintf(1, 'Reading file number: %d   name: %s\n', (fileNum - startFileNum + 1), fileName);
    ADCData = fread(fid, 'int16=>int16');
    fclose(fid);

    RMSRecordLength = 60.0 * Fs;    % samples
    RMSIndex = 1;
    processRMS = 1;
    if( processRMS == 1)
        RMSRecordNum = 1;
        startSample = ((RMSRecordNum - 1) * RMSRecordLength) + 1;
        lastSample = startSample + RMSRecordLength - 1;
        while( lastSample <= samplesPerFile)
            RMS(RMSIndex) = sqrt(mean(single(ADCData(startSample:lastSample)).^2));
            fprintf('RMS index %d  RMS Value: %f\n', RMSIndex, RMS(RMSIndex));
            RMSRecordNum = RMSRecordNum + 1;
            RMSIndex = RMSIndex + 1;
            startSample = ((RMSRecordNum - 1) * RMSRecordLength) + 1;
            lastSample = startSample + RMSRecordLength - 1;
        end
    end
    
    histRecordLength = 2 * 60.0 * Fs;   % samples
    numHistBins = 21;
    binSpacing = (2^16) / numHistBins;
    histX = (-(numHistBins - 1) *  binSpacing / 2):binSpacing:((numHistBins - 1) * binSpacing ) / 2;
    histIndex = 1;
    processHist = 1;
    if( processHist == 1)
        histRecordNum = 1;
        startSample = ((histRecordNum - 1) * histRecordLength) + 1;
        lastSample = startSample + histRecordLength - 1;
        while( lastSample <= samplesPerFile)
            histY(histIndex,:) = hist(single(ADCData(startSample:lastSample)), histX);
            fprintf('Process Histogram index %d\n', histIndex);
            histRecordNum = histRecordNum + 1;
            histIndex = histIndex + 1;
            startSample = ((histRecordNum - 1) * histRecordLength) + 1;
            lastSample = startSample + histRecordLength - 1;
        end
    end

    PSDSegmentLength = 2^16;    % samples
    PSDSegmentsPerRecord = 8;
    PSDRecordLength = PSDSegmentsPerRecord * PSDSegmentLength;
    overlapPercent = 50;
    hSpectrum = spectrum.welch('Blackman-Harris', PSDSegmentLength, overlapPercent);
    PSDIndex = 1;
    processPSD = 0;
    if( processPSD == 0)
        PSDRecordNum = 1;
        startSample = ((PSDRecordNum - 1) * PSDRecordLength) + 1;
        lastSample = startSample + PSDRecordLength - 1;
        VData = single(LSB * single(ADCData));
        while( lastSample <= samplesPerFile)
            pxx = psd(hSpectrum, single(VData(startSample:lastSample)), 'Fs', Fs);
            Pxx(PSDIndex,:) = 10 * log10(pxx.Data); %turn the linear pxx values into logarithmic dB values
            fprintf('Process PSD index %d\n', PSDIndex);
            PSDRecordNum = PSDRecordNum + 1;
            PSDIndex = PSDIndex + 1;
            startSample = ((PSDRecordNum - 1) * PSDRecordLength) + 1;
            lastSample = startSample + PSDRecordLength - 1;
        end
    end
    clear ADCData VData
end

procDataFileName = strcat('procData',datestr(now),'.mat');
procDataFileName = strrep(procDataFileName, '-', '');
procDataFileName = strrep(procDataFileName, ':', '');
procDataFileName = strrep(procDataFileName, ' ', '-');
fprintf(1, 'Writing Matlab data file %s\n', procDataFileName);
save(procDataFileName);

figure
plot(RMS);
figure
waterfall(histX,1:histIndex-1,histY)


tElapsed = toc(tStart)


