% Frist run 02_21_10 - Bill Ussler
% 5 minute data interval between files

clear all
close all

%Input actual files size used on deployment
samplesPerFile = 20.0 * 60 * 48000.0;   %5 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\'; % windows - format only
%dirName = '/Volumes/Larry/Deployment_Jan_22_10/R2010-021-16-28-19/R00000/'; % Bill Mac
startFileNum = 4300;
numFilesToProcess = 50;

endFileNum = startFileNum + numFilesToProcess - 1;

tStart = tic;

RMSIndex = 1;
histIndex = 1;
PSDIndex = 1;

histIndex = 1;
for fileNum=startFileNum:endFileNum
    fileName = strcat(dirName, 'D',  sprintf('%06d', fileNum));
    fid = fopen(fileName);
    fprintf(1, '\nReading file number: %d   name: %s\n', (fileNum - startFileNum + 1), fileName);
    ADCData = fread(fid, 'int16=>int16');
    fclose(fid);
    
    RMSRecordTime = 60.0; % seconds
    RMSRecordLength = RMSRecordTime * Fs;    % samples
    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
        RMSTimeIndex = (0:RMSIndex-2)*RMSRecordTime; % time index in seconds      
    end
    
    histRecordTime = 60.0; % seconds
    histRecordLength = histRecordTime * 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
        histTimeIndex = (0:histIndex-2)*histRecordTime; % time index in seconds
    end

    PSDSegmentLength = 2^16;    % samples
    PSDSegmentsPerRecord = 64;
    PSDRecordLength = PSDSegmentsPerRecord * PSDSegmentLength;
    PSDRecordTime = PSDRecordLength/Fs; 
    overlapPercent = 50;
    hSpectrum = spectrum.welch('Blackman-Harris', PSDSegmentLength, overlapPercent);
    
    processPSD = 0;
    if( processPSD == 1)
        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
        PSDTimeIndex = (0:PSDIndex-2)*PSDRecordTime; % time index in seconds
    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(1)
plot(RMSTimeIndex,RMS);
xlabel('Time (sec)')
ylabel('RMS')
title('January 22, 2010 Deployment')
figure(2)
waterfall(histX,histTimeIndex,histY)
xlabel('Amplitude')
ylabel('Time (sec)')
zlabel('Number of Samples')
title('January 22, 2010 Deployment')
%figure(3)
%waterfall(pxx.frequencies,PSDTimeIndex,Pxx)
%xlabel('Frequency (Hz)')
%ylabel('Time (sec)')
%zlabel('Amplitude')

tElapsed = toc(tStart)


