% tests on data transfer for channel sounder

format
format compact
clear

blocks = input('Enter number of 1K blocks per dataset (16 or 128): ');


dataset_size = 1024 * blocks;      % this is the size of the dataset stored on disc
source_increment = 1;      % this is the distance between samples on the disc
plot_size = 1024 * blocks;         % this is the size we're going to plot

% This program plots FFTs from the data collected on the TDS PCMCIA Card
% It is only a 'development' version and is not 'user friendly'.
% Please set above parameters as required.

word_size = 4;
skip_distance = dataset_size / plot_size - 1;

% This is the distance (-1) between the samples we use from the disc.
% Now calculate the axis labelling

effective_sample_frequency = 500000/(source_increment * (skip_distance+1));
resolution = effective_sample_frequency/plot_size;
freq_axis = effective_sample_frequency/1000;

% Now open the disc file and skip over the initial (previously logged) data.

skip_bytes = skip_distance * word_size
fid3=fopen('c:\my documents\chs_data\data.fil','r');

skipover = input('Enter hex number for skipover distance  $','s');
skipover_bytes = hex2dec(skipover)-4;
skip = fread(fid3,1,'integer*4',skipover_bytes);

% Now read a dataset into f3 and skip the 'post-amble'

f3 = fread(fid3,plot_size,'integer*4',skip_bytes);
skip = fread(fid3,2048,'integer*4',skip_bytes);
first_word_f3 = checkhex(f3(1))
averaging3 = bitand(hex2dec(checkhex(f3(1))),255);

% Now read a dataset into f4, in similar fashion

f4 = fread(fid3,plot_size,'integer*4',skip_bytes);
skip = fread(fid3,2048,'integer*4',skip_bytes);
first_word_f4 = checkhex(f4(1))
averaging4 = bitand(hex2dec(checkhex(f4(1))),255);

fclose('all');

% 'averaging3' holds the signal-averaging setting. 
% We process the data as follows...
%
% - remove DC offset from signal, and scale to +/- 1.25V
% - divide signal by signal-averaging setting
% - calculate mean and rms
% - calculate DFT

[f3, F3, meanf3, rmsf3] = process_dataset(f3,averaging3);

[f4, F4, meanf4, rmsf4] = process_dataset(f4,averaging4);

meanf3
rmsf3
meanf4
rmsf4

% ---------------------------------------------------

% Plot f3 samples in time domain and frequency domain

figure(1)
plot(f3);
title(strcat('time-domain:', ...
   ' Mean(pre-zero)=',num2str(meanf3),'mV. RMS=',num2str(rmsf3),'mV. sets=',num2str(averaging3),'.'))

Vresolution = 1E6 * 1.25 / 2^15 / averaging3;

figure(2)
plot(linspace(0,freq_axis,length(F3)),abs(F3))
title(strcat('freq-domain (kHz):', ...
   ' Resolution=',num2str(resolution),'Hz/',num2str(Vresolution),'uV. sets=',num2str(averaging3),'.'))

% NB 38uV is 1.25V / 2^15

ax = axis;
axis([0 freq_axis/2 ax(3:4)])

% ---------------------------------------------------

% Plot f4 samples in time domain and frequency domain

figure(3)
plot(f4);
title(strcat('time-domain:', ...
   ' Mean(pre-zero)=',num2str(meanf4),'mV. RMS=',num2str(rmsf4),'mV. sets=',num2str(averaging4),'.'))

Vresolution = 1E6 * 1.25 / 2^15 / averaging4;

figure(4)
plot(linspace(0,freq_axis,length(F4)),abs(F4))
title(strcat('freq-domain (kHz):', ...
   ' Resolution=',num2str(resolution),'Hz/',num2str(Vresolution),'uV. sets=',num2str(averaging4),'.'))

axis([0 freq_axis/2 ax(3:4)])

% ---------------------------------------------------

% Zoom in, frequency domain

figure(5)
plot(linspace(0,freq_axis,length(F3)),abs(F3))
axis([190 200 ax(3:4)])

Vresolution = 1E6 * 1.25 / 2^15 / averaging3;
title(strcat('freq-domain (kHz):', ...
   ' Resolution=',num2str(resolution),'Hz/',num2str(Vresolution),'uV. sets=',num2str(averaging3),'.'))

% ----------------------------------------

% to check the calibration, 
% we replace f3 with a square wave of amplitude +/- 2^23.
% This ought to display as +/- 1250mV in the time domain
% and as 4/pi times this in the freq domain, i.e. +/- 1590mV

cal(1:32)=2^23; 
cal(33:64)=-(2^23);
for j=6:13
   cal(2^j+1:2^(j+1)) = cal(1:2^j);
end

fcal = reshape(cal,length(cal),1);

[fcal, Fcal, meanfcal, rmsfcal] = process_dataset(fcal,1);

meanfcal
rmsfcal

figure(6)
plot(fcal);
title(strcat('time-domain:', ...
   ' Mean(pre-zero)=',num2str(meanfcal),'mV. RMS=',num2str(rmsfcal),'mV. sets=',num2str(1),'.'))

Vresolution = 1E6 * 1.25 / 2^15 ;

figure(7)
plot(linspace(0,freq_axis,length(Fcal)),abs(Fcal))
title(strcat('freq-domain (kHz):', ...
   ' Resolution=',num2str(resolution),'Hz/',num2str(Vresolution),'uV. sets=',num2str(1),'.'))

% NB 38uV is 1.25V / 2^15

ax = axis;
axis([0 freq_axis/2 ax(3:4)])

% ---------------------------------------------------

% to display as hex we need to convert to unsigned

%h = 2^32 - 1 - f;
%h = bitcmp(bitand(h,2^32-1),32);
%dec2hex(h(1),8)  % should be FFAA0001
%dec2hex(h(2),8)  % should be 00275A01
%dec2hex(h(16384),8)  % should be 00795401

% pick out the peaks
%for i=1:length(F3),
%   if (abs(F3(i)) > 100)
%      disp([abs(F3(i)),i])
%   end
%end
