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

clear all
close all

runStartTime = datestr(now);
rootProcDataFileName = strcat('procData', runStartTime);
rootProcDataFileName = strrep(rootProcDataFileName, '-', '');
rootProcDataFileName = strrep(rootProcDataFileName, ':', '');
rootProcDataFileName = strrep(rootProcDataFileName, ' ', '-');
rootProcDataFileName = strcat(rootProcDataFileName, '-');
procDataFileNum = 1;

%Input actual files size used on deployment
samplesPerFile = 5.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 = 'Q:\DeploymentJan22-2010\R2010-021-16-28-19\R00000\'; % windows - format only
%dirName = '/Volumes/Larry/Deployment_Jan_22_10/R2010-021-16-28-19/R00000/'; % Bill Mac
startFileNum = 258;
numFilesToProcess = 10;

endFileNum = startFileNum + numFilesToProcess - 1;

tStart = tic;

RMSIndex = 1;
histIndex = 1;
PSDIndex = 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);
    if(fileNum == startFileNum)
        fileData = dir(fileName);
        startDatenum = fileData.datenum + (1/3);    % add 8 hours to the time stamp to correct for OS misinterpretation
    end
    
    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
     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;
    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 = 64;
    PSDRecordLength = PSDSegmentsPerRecord * PSDSegmentLength;
    PSDRecordTime = PSDRecordLength/Fs; 
    overlapPercent = 50;
    hSpectrum = spectrum.welch('Blackman-Harris', PSDSegmentLength, overlapPercent);
    processPSD = 1;
    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
    end
    clear ADCData VData
    
    if( mod( ((fileNum - startFileNum) + 1), 864) == 0 )
        procDataFileName = strcat(rootProcDataFileName, num2str(procDataFileNum), '.mat');
        fprintf(1, 'Writing Matlab data file %s\n', procDataFileName);
        save(procDataFileName);
        clear Pxx 
        procDataFileNum = procDataFileNum + 1;
    end
end

if(exist('Pxx'))
    procDataFileName = strcat(rootProcDataFileName, num2str(procDataFileNum), '.mat');
    fprintf(1, 'Writing Matlab data file %s\n', procDataFileName);
    save(procDataFileName);
end

% if(processRMS == 1)
%     RMSTimeIndex = startDatenum + (((0:RMSIndex-2)*RMSRecordTime) / 86400.0); % time index in seconds
%     figure
%     plot(RMSTimeIndex,RMS);
%     datetick('x')
%     xlabel('Clock Time')
%     ylabel('RMS Amplitude')
%     title('January 22, 2010 Deployment')
% end
% 
% if(processHist == 1)
%     figure
%     histTimeIndex = startDatenum + (((0:histIndex-2)*histRecordTime) / 86400.0); % time index in serial date number format
%     waterfall(histX,histTimeIndex,histY)
%     datetick('y')
%     xlabel('ADC Counts')
%     ylabel('Clock Time')
%     zlabel('Number of Samples')
%     title('January 22, 2010 Deployment')
% end
% 
% if(processPSD == 1)
%     figure
%     PSDTimeIndex = startDatenum + (((0:PSDIndex-2)*PSDRecordTime) / 86400.0); % time index in serial date number format
%     surf(pxx.frequencies, PSDTimeIndex, double(Pxx), 'EdgeColor', 'none')
%     datetick('y')
%     xlabel('Frequency (Hz)')
%     ylabel('Clock Time')
%     zlabel('Amplitude')
% end
        
tElapsed = toc(tStart)


