
% Determin the optimal ratios for 2011 Antarctica TO from recorded waveforms in the field
% Chebshev window is applied in frequency domain
% Calculate the PSL and ISL
% Jilu Li,Nov. 25,2011

clear all; 
close all;
clc;

%1) Reference chirp signal 
f0=69e6;      % Linear chirp starting frequency,Hz
f1=99e6;      % Linear chirp ending frequency,Hz
fc=(f0+f1)/2;
BW=f1-f0;      % Linear chirp bandwidth,Hz  
LO=84e6;
pulseDuration=1e-6;          %linear chirp duration,sec
alpha=pi*BW/pulseDuration;    %Chirp rate
wc=2*pi*fc;
w0=2*pi*f0;
Fs=111e6;     % Sampling Frequency
dt=1/Fs;
[B A] = butter(8,BW/Fs*1.1);
n = pulseDuration/dt + 1;           % number of time samples
t=(0:n-1)'*dt;                      % time array for data acquisition
s=cos(w0*t+alpha*(t.^2));           % waveform
sw = s.*chebwin(size(s,1),100);     % windowed waveform
taperR =0:0.01:0.3;
PSL = zeros(31,1);
ISL = zeros(31,1);
pulseWidth = zeros(31,1);
for ii = 1:31
staper = s.*tukeywin(n,taperR(ii));        % amplitude tapering

%2) Baseband conversion
sb = s.*exp(-j*2*pi*LO*t);
sb = filtfilt(B,A,sb);
sbf= fftshift(fft(sb));
staperb = staper.*exp(-j*2*pi*LO*t);
staperb = filtfilt(B,A,staperb);
staperbf = fftshift(fft(staperb));
df=1/t(end);
f=(-n/2:n/2-1)*df;
ifirst=find(f>=-BW/2,1,'first');
ilast=find(f<=BW/2,1,'last');
sbfw = sbf;
sbfw(ifirst:ilast)=sbf(ifirst:ilast).*chebwin(ilast-ifirst+1);
sbfw(1:ifirst-1)=0;
sbfw(ilast+1:end)=0;
staperbfw = staperbf;
staperbfw(ifirst:ilast)=staperbf(ifirst:ilast).*chebwin(ilast-ifirst+1);
staperbfw(1:ifirst-1)=0;
staperbfw(ilast+1:end)=0;

%3)match filtering 
sm = fftshift(ifft(sbf.*conj(sbfw)));
smNormal = 20*log10(abs(sm/max(abs(sm))));
smtaper = fftshift(ifft(staperbf.*conj(staperbfw)));
smtaperNormal =20*log10(abs(smtaper/max(abs(smtaper))));
% figure;plot(smtaperNormal);
% pause

%4) calculate pulse width,PSL and ISL
[tmp,mainPeakInd] = max(smNormal);
firstNullInd3dB = find(smNormal(mainPeakInd:end)<=-3,1,'first');
pulseWidth0 = 2*(firstNullInd3dB-1)*dt*BW;
[tmp,mainPeakInd] = max(smtaperNormal);
smtaper_diff = diff(smtaperNormal(mainPeakInd:end));
firstNullInd3dB = find(smtaperNormal(mainPeakInd:end)<=-3,1,'first');
firstNullInd = find(smtaper_diff>=0,1,'first');
pulseWidth(ii) = 2*(firstNullInd3dB-1)*dt*BW/pulseWidth0;
PSL(ii) = max(smtaperNormal(mainPeakInd+firstNullInd:end));
ISL(ii) = 10*log10(sum(abs(smtaper(mainPeakInd+firstNullInd:end)).^2)/(sum(abs(max(smtaper(mainPeakInd:mainPeakInd+firstNullInd-1))))^2));
end

figure;plot(taperR,PSL)
set(gca,'fontSize',18)
axis([0.05,0.3,-100 0]);grid;
hold on;plot(taperR,ISL,'r--')
legend('PSL','ISL')
xlabel('taper ratio')
ylabel('PSL,ISL(dB)')
% print -dpdf PSLISL10us