%
% resampANDfilt_example.m
%
% Example file presenting the resample function, PSD calculation
% and filtering.
%

% Set up chanvector for filelist operation
[gpsTimes,frameFiles,frameDurs]=dir2framelist('./samrds/S3/L3/LHO/H-RDS_R_L3-7545');

% Get data using chanvector
channel = chanstruct('H1:LSC-AS_Q');
[data, fs] = chanvector(channel, 754562512, 300,gpsTimes,frameFiles,frameDurs); %300 seconds of S3 playground data

% Downsample from 16384Hz to 4096Hz
% make a plot showing before and after

[psdBefResamp, fBefResamp] = psd(data, 25600, fs, [], 25600/2);
% psd input and output arguments:
%    psdBefResamp - psd
%    fBefResamp   - frequencies corresponding to each value of the psd vector
%    data         - data that you are calculating the psd for
%    25600        - number of frequencies that the psd is calculated at 
%                   (length of psdBefResamp and fBefResamp)
%    fs           - sample rate of data
%    []           - window used in psd estimation, [] gives the defult window
%                   (hanning window of the same length as the psd)
%    25600/2      - length of overlap

figure
subplot(2,1,1)
semilogx(fBefResamp, 20*log10(psdBefResamp))
titleString = sprintf('%s PSD of H1:LSC-AS\\_Q Before and After Resampling', datestr(now));
title(titleString);
ylabel('PSD (dB)');
legend('Before Resampling');
axis([2 16384 -80 100])

% Downsample to 4096Hz
% The zero padding removes the delay caused by resample
initialFreq = fs;
finalFreq = 4096;
[p,q] = rat(finalFreq/initialFreq); % ratio, in this case p = 1, q = 4
padData = [0 data];
resampData = resample(padData,p,q);
resampData = resampData(2:end);
resampFs = initialFreq*p/q;

[psdAftResamp, fAftResamp] = psd(resampData, 2560, resampFs, [], 2560/2);
subplot(2,1,2);
semilogx(fAftResamp, 20*log10(psdAftResamp))
ylabel('PSD (dB)');
legend('After Resampling');
xlabel('Frequency (Hz)')
axis([2 16384 -100 100])

% Filter to suppress low frequencies
% Make a plot showing data before and after filtering

% Input arguments to remez are as follows:
% remezFilter = remez(filtLen, edges, amps, weights);
% filtLen = length of remezFilter
% edges = transition points of filter response as a fraction of Nyquist
% amps = desired amplitudes of filter response at edges
% weights = how much the error is minimized between pairs of edges relative 
%   to everywhere else
% these inmput arguments found via remezord
filterCoeffs = remez(138, [0 0.0244 0.0625 1]', [0 0 1 1]', [10 1]');
% this filter supresses frequencies below 128 Hz
figure
filteredData = filter(filterCoeffs,1,resampData);
[psdAftFilt, fAftFilt] = psd(filteredData, 2560, resampFs, [], 2560/2);
plot(fAftResamp, 20*log10(psdAftResamp), fAftFilt, 20*log10(psdAftFilt))
titleString = sprintf('%s PSD of H1:LSC-AS\\_Q Before and After Filtering', datestr(now));
title(titleString);
ylabel('PSD (dB)');
xlabel('Frequency (Hz)');
axis([0 200 -180 100])
legend('Before Filtering', 'After Filtering');