% blimit_example.m
%
% Example file displaying base-banding
%
% Contains hard coded band parameters

band.freqs=[384 512]; % Band to move to zero frequency
%  low-pass filter for this band
band.blp = remez(32, [0 .5 .6250 1]', [1 1 0 0]', [1 1]'); 
                                                           
  % 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');

% 300 seconds of S3 playground data
[data, fs] = chanvector(channel, 754562512, 300,gpsTimes,frameFiles,frameDurs); 

% Produce Psd before
[psdBefBband, fBefBband] = psd(data, 51200, fs, [], 51200/2);

% Find range of indices spaning the frequency band
fRangeBef=find( fBefBband>=band.freqs(1) & fBefBband<=band.freqs(2) );

figure
subplot(2,1,1)
plot(fBefBband(fRangeBef), 20*log10(psdBefBband(fRangeBef)))
titleString = sprintf('%s PSD of H1:LSC-AS\\_Q Before and After Base-banding', datestr(now));
title(titleString);
ylabel('PSD (dB)');
xlabel('Frequency (Hz)');
legend('Before Base-banding');
xlim([fBefBband(fRangeBef(1)) fBefBband(fRangeBef(end))]);

% Baseband using blimit (uses fir1 and resample)
[dataBband,fsBband] = blimit(data,fs,band);

[psdAftBBand, fAftBband] = psd(dataBband, 512, fsBband, [], 512/2);
fRangeAft=find( fAftBband>=0 & fAftBband<=(band.freqs(2)-band.freqs(1)) );

subplot(2,1,2);
plot(fAftBband(fRangeAft), 20*log10(psdAftBBand(fRangeAft)))
ylabel('PSD (dB)');
xlabel('Frequency (Hz)');
legend('After Base-banding');
xlim([fAftBband(fRangeAft(1)) fAftBband(fRangeAft(end))]);