% investigation of power capability when using sequence 
% with untuned antenna

filetype = '-dtiff';

N = 256;
f = [ ...
  1 0 1 1 0 0 0 1 1 1 1 0 1 0 0 0 0 1 1 1 1 1 1 1 1 0 0 1 0 0 0 0 ...
  1 0 0 ...
        1 0 0 1 1 1 1 1 0 1 0 1 0 1 0 1 1 1 0 0 0 0 0 1 1 0 0 0 1 ...
  0 1 0 1 1 0 0 1 1 0 0 1 0 1 1 1 1 1 1 0 1 1 1 1 0 0 1 1 0 1 1 1 ...
  0 1 1 1 0 0 1 0 1 0 1 0 0 1 0 1 0 0 0 1 0 0 1 0 1 1 0 1 0 0 0 1 ...
  1 0 0 1 1 1 0 0 1 1 1 1 0 0 0 1 1 0 1 1 0 0 0 0 1 0 0 0 1 0 1 1 ...
  1 0 1 0 1 1 1 1 0 1 1 0 1 1 1 1 1 0 0 0 0 1 1 0 1 0 0 1 1 0 1 0 ...
  1 1 0 1 1 0 1 0 1 0 0 0 0 0 1 0 0 1 1 1 0 1 1 0 0 1 0 0 1 0 0 1 ...
  1 0 0 0 0 0 0 1 1 1 0 1 0 0 1 0 0 0 1 1 1 0 0 0 1 0 0 0 0 0 0 0 ];

% ---------------------------------  alternatively use M-seq
%f = [ f(1:34) f(36:256) ];  N = 255;

% ---------------------------------  alternatively use square wave
while 0,
f = [0 1 0 1]; f = [f f f f]; f = [f f f f]; f = [f f f f];
N = length(f);
% this results in current of +/- 0.5 and power of 1/12
end

% Expand the time axis, as this gives us a crude method of calculating
% the power, without resorting to formulas as we did in the old BASIC version
% of this program.

% First make f have 'scale' identical rows
scale = 64;
f = [f;f;f;f];
f = [f;f;f;f];
f = [f;f;f;f];

% now reshape it into a single row
N = scale*N;
f = reshape(f,1,N);

% Calculate integral of f, which is current waveform

current(1)=0;
for k=2:N,
    current(k) = current(k-1) + (2*f(k-1) - 1)/scale;
end

% make sure current vector is the right length
current = current(1:N);

imax = max(current);
imin = min(current);

% calculate mean current, and bias the current for zero mean
current = current - sum(current)/N;
%current = current - (imax+imin)/2;

% shift current so that plot begins at zero
cut = 86*64+1;
current = [current(cut:length(current)) current(1:cut-1)];

% print current relative to square wave value
imax = max(current)/0.5
imin = min(current)/0.5
irange = imax-imin

%calculate and print mean power relative to sq wave value
power = sum(current.^2)/N / (1/12)


% -------------------------------- plot integrated sequence ----
fig=1; figure(fig); clf
plot(linspace(1,256,16384),current, 'b-','linewidth',2);
set(gca,'Xtick',linspace(32,256,8))
hold on; plot([1 256],[0 0])
axis tight
set(gca,'LineWidth',1.5,'FontSize',12,'TickDir','out')
xlabel('Sample Number','fontsize',12)
%ylabel('Integrated Sequence','fontsize',12)
title('Integrated Channel-Sounding Sequence','fontsize',12)
print(filetype, '', ['seq-' num2str(fig)] )

% -------------------------------- plot spectrum ----
fig=2; figure(fig); clf
F = fft(2*f-1)/8192;
semilogy(abs(F).^2, 'b-','linewidth',2);
set(gca,'Xtick',linspace(0,4*256,5))
axis([0 4*256 0.000001 0.2])
set(gca,'LineWidth',1.5,'FontSize',12,'TickDir','out')
xlabel('Sample Number','fontsize',12)
ylabel('Power Spectrum','fontsize',12)
title('Fourier Transform of Sequence','fontsize',12)
print(filetype, '', ['seq-' num2str(fig)] )


