echo off
clc
%	PREDISTORTION INTERMODULATION STUDY    predis1.m
%
% Simple case of amplitude distortion only.
%
% RSL  15 Nov 00

axis([1 2 3 4]);  axis;			% Auto scale plots

n = 512;
ndb = 20*log10(n/2);
pi2 = 2.0*pi;

% Amplifier:
c1=1.0;
c3=-0.1;
% c5=-0.05;
c5=0.0;
% Assuming a 1 Ohm system
ip3=10.0*log10(-0.66666667*c1*c1*c1/c3);
fprintf('\nPredistortion IM Polynomial predistorter\n');
fprintf('All signals are in a 1 Ohm system,');
fprintf(' and power units are Watts or dBW\n\n');
fprintf('c1=%.3f   c3=%.3f   IP3=%.3f dBW', c1, c3, ip3);
fprintf('   c5=%.3f\n\n', c5);

% Input signals:
a1 = input('Max sine wave (each one) amplitude A1, as in A1*sin(x)? ');
kf1=pi2*17.0;
a2 = a1;
kf2=pi2*23.0;

% Power fundamental, w/o compression
pp1=10.0*log10(0.5*a1*a1);
% From power series:
ppim3=3*pp1-2*ip3;
im3_suppr_dB=pp1-ppim3;
fprintf('Each signal if there was no compression, Pwr=%.3f dBW\n', pp1);
fprintf('Cubic power series Amp IM3 supr=%.3f dB\n\n', im3_suppr_dB);
pause

% Sample rate 512 per sec:
t = 0:0.001953125:1.171875;		%  <-- 0 to 600 data points

% For two signals the phase should be irrelevant:
ph1=pi2*0.22;
ph2=pi2*0.87;

si=a1*cos(ph1+kf1*t);
si = si + a2*cos(ph2+kf2*t);

% Plot time waveform:
axis([0 0.2 -2.0 2.0]);
plot(t(1:256),si(1:256)), title('Waveform, si'), ..
xlabel('time'),grid
pause;
clg;

% Predistort:
for k=1:1:512
   % Absolute value 5th order predistorter:
%  a=abs(si(k));
%  g= 0.9999+0.0018*a+.0891*a*a+.0251*a*a*a-0.0155*a*a*a*a+0.0107*a*a*a*a*a;
   %
   % Quartic predistorter:
   b=si(k);
   g=1.0147-0.0409*b*b+0.1930*b*b*b*b;
%   si(k)=g*si(k);
end;

% Plot time waveform:
axis([0 0.2 -2.0 2.0]);
plot(t(1:256),si(1:256)), title('After predistortion, si'), ..
xlabel('time'),grid
pause;
clg;


% Amplify, si, so are amplifier input and output:
for k=1:1:512
    so(k) = c1*si(k) + c3*si(k)*si(k)*si(k);
    so(k) = so(k) + c5*si(k)*si(k)*si(k)*si(k)*si(k);
end;

% Plot time waveform:
axis([0 0.2 -2.0 2.0]);
plot(t(1:256),so(1:256)), title('Amplifier output, so'), ..
xlabel('time'),grid
pause;
clg;


si=[];					% Reclaim memory

Y = fft(so(1:512), 512);		% Output spectra of amp
Pyy = -3.01+10*log10(Y.*conj(Y));	% -3.01 puts into dBW
Prel = Pyy-Pyy(18)-6.02;		% -6.02 for peak thing

fprintf('DFT based Amp IM3 supr=%.3f dB\n\n', Pyy(18)-Pyy(12));

% Note that  spectrum for f=17 Hz is Pyy(18), etc.

%	To plot the power spectrum, we must first form a frequency axis:
f = (0:512);
axis([0 50 -60 0]);
plot(f(1:128),Prel(1:128)), title('Amplifier Relative Output Spectrum'), ..
xlabel('Frequency (Hz)'),grid
pause;
clg;
for k=16:21
    if Pyy(k) <=-100
	Pyy(k)=-99.9;
    end;
end;
fprintf('\nAmp dBW out: f=5, 11, 17, 23, 29, 35: %.1f %.1f',Pyy(6),Pyy(12) );
fprintf(' %.1f %.1f',Pyy(18),Pyy(24) );
fprintf(' %.1f %.1f',Pyy(30),Pyy(36) );