%Run from editor debug(F5)-DPSK UWB system analysis
%JC 12/10/05
%This m file simulates a differential phase shift keyed (DPSK) ultra wide bandwidth(UWB) system using
%a fifth derivative waveform equation of a Gaussian pulse. I thought the setup for DPSK
%would be apparent to configure from my other fifth drivative m files but several emails 
%convinced me otherwise so here it is. This system should be fairly easy to breadboard in
%the lab assuming an appropriately starting pulse can be realized along with the required 
%delay line in the receuver. I have to keep reminding myself we are only talking one to
%five meters, not miles. One could walk in front of the antennas and flip and move the antennas to
%observe multipath, etc. Another interesting problem for test equipment vendors
%is how to measure and characterize the spectrum of the waveform coming off the 
%transmitting antenna in the far field. Anyway, good luck to those who like to tinker and
%play in new areas of wireless. 
%The receiver section uses a squaring %circuit(mixer) with a delay in one leg and a low pass filter
%is used to integrate the waveforms to recover the voltage levels(a radiometer and/or correlaton receiver).
%A one bit flash ADC is simulated via a comparator routine to show the demodulated 
%output bits from the integrator. Various plots of time and frequency waveforms are shown at
%different positions in the file so the user can observe and analyze the
%operation of the simulated system. The transmitter is simulated assuming
%the fifth order waveforms are being sent from an appropriate antenna. The receiving antenna
%in a real system would distort the original and delayed signal equally. 
%The fifth derivative waveform equation is obtained by use of the Matlab symbolic processor in
%matlab.You would actually use the fourth derivative of the gaussian
%monocycle 1st derivative which is t/pw*exp(-t^2/pw^2).Procedure as follows:
%syms t pw 
%t/pw*exp(-t^2/pw^2)enter
%diff(ans,t,4)
%ans=60/pw^5*t*exp(-t^2/pw^2)-80/pw^7*t^3*exp(-t^2/pw^2)+16*t^5/pw^9*exp(-t^2/pw^2)
%set pw=1 to normalize and do an ezplot(ans) and you should get a plot of the fifth order
%derivative which has 5 zero crossings.It is apparent from my previous files that the 1st and
%second derivative pulses will not meet the FCC spectral mask without
%reducing the transmitter power. The fifth order derivative pulse meets the FCC
%mask without a reduction in power. It also retains a wide 3DB and 10DB bandwidth.
%One may ask how to generate these pulses at the transmitter. One possible
%method might be to take a high speed differentially encoded digital bit stream and process it
%thru highpass or bandpass filters(the pulse would ring on the leading edge
%and would ring and be inverted on the trailing edge) with the output being a close match to a 
%fourth derivative waveform and use a TX antenna that differentiates the pulse
%once for fifth order. DPSK does not require a template at the receiver and any 
%distortions caused by the receiving antenna or the channel path would be reduced
%and possibly balanced out. It may be possible to improve the system (for short distances 
%and very high bit rates), to use direct sequence spread spectrum(DSSS) (for smoothing and 
%multiple access-may not be required since UWB is already a spread spectrum process) and forward 
%error correction(FEC) coding with the DPSK scheme. Not needing to synch up the receiver with a templet
%solves a lot of design problems. You would only lose several DB of Eb/No(1e-3) 
%and as stated before would not require generation of a template at 
%the receiver for syncronization assuming timing jitter can be held to a minimum.
%The characteristics of the fifth order waveform are as follows and which you can verify are:
%pw(tail to tail)=~.5e-9
%fc=~7Ghz
%3DB fl to fh=~3GHz
%10DB fl to fh=~6GHz
%You might wish to review my other files published under UWB or author(C,J) to get
%additional information on UWB analysis.
%================================================
clear
Fs=100e9;%sample frequency
Fn=Fs/2;%Nyquist frequency
t=-1e-9:1/Fs:46e-9;%time vector sampled at Fs Hertz. zoom in/out using (-1e-9:1/Fs:xxxx)-set for
%6e-9 to see first pulse.
%================================================ 
% EQUATIONS
%================================================
%y=A*(t/pw).*exp(-(t/pw).^2);%1st derivative of Gaussian pulse=Gaussian monocycle
%y =1*(1 - 4*pi.*((t)/pw).^2).* exp(-2*pi.*((t)/pw).^2);%2nd derivative of Gaussian
%pulse=doublet(two zero crossings)
pw=75e-12;%value sets pulse width of fifth derivative looking at figure 1
y=((60./pw.^5).*(t-0).*exp(-(t-0).^2./pw.^2)-(80./pw.^7).*(t-0).^3.*exp(-(t-0).^2./pw.^2)+ ...
(16.*(t-0).^5./pw.^9.*exp(-(t-0).^2/pw.^2)));%Fifth derivative of Gaussian pulse.(5 zero crossings)

yd=((60./pw.^5).*(t-5e-9).*exp(-(t-5e-9).^2./pw.^2)-(80./pw.^7).*(t-5e-9).^3.*exp(-(t-5e-9).^2./pw.^2)+ ...
(16.*(t-5e-9).^5./pw.^9.*exp(-(t-5e-9).^2/pw.^2)));%(one) fifth derivative pulse delayed by 5e-9 sec
%================================================
%NOISE SETUP FOR BER AND SNR
%================================================

noise=1e-50*(randn(size(t)));%(Noise-AWGN)Set to 1e-50 to disable(set to~3e-5 for Eb/No=9DB BER @1e-3) 

%================================================
%DPSK  
%================================================
%The following series of equations sets the pulse recurring frequency(PRF)
%at 200MHz(waveform repeats every 5e-9 sec and a
%modulated bit stream(bit rate=200Mb/s)of 1001101 (differentially encoded to a value of 01110110) 
%where a phase change (0to1 or 1to0)=0 and no phase change (0to0 or 1to1)=1
 
%==================================================
% FIFTH DERIVATIVE(DPSK) WITH 5 PULSES)
%==================================================

%DPSK modulated fifth(ydiff)
A=.6e-45;%sets voltage level out of TX or input to mixer(.6e-45 for .6mv volt peak to peak)
ydiff=-A*y+ ...
A*((60./pw.^5).*(t-5e-9).*exp(-(t-5e-9).^2./pw.^2)-(80./pw.^7).*(t-5e-9).^3.*exp(-(t-5e-9).^2./pw.^2)+ ...
(16.*(t-5e-9).^5./pw.^9.*exp(-(t-5e-9).^2/pw.^2)))+ ...
A*((60./pw.^5).*(t-10e-9).*exp(-(t-10e-9).^2./pw.^2)-(80./pw.^7).*(t-10e-9).^3.*exp(-(t-10e-9).^2./pw.^2)+ ...
(16.*(t-10e-9).^5./pw.^9.*exp(-(t-10e-9).^2/pw.^2)))+ ...
A*((60./pw.^5).*(t-15e-9).*exp(-(t-15e-9).^2./pw.^2)-(80./pw.^7).*(t-15e-9).^3.*exp(-(t-15e-9).^2./pw.^2)+ ...
(16.*(t-15e-9).^5./pw.^9.*exp(-(t-15e-9).^2/pw.^2)))+ ...
-A*((60./pw.^5).*(t-20e-9).*exp(-(t-20e-9).^2./pw.^2)-(80./pw.^7).*(t-20e-9).^3.*exp(-(t-20e-9).^2./pw.^2)+ ...
(16.*(t-20e-9).^5./pw.^9.*exp(-(t-20e-9).^2/pw.^2)))+ ...
A*((60./pw.^5).*(t-25e-9).*exp(-(t-25e-9).^2./pw.^2)-(80./pw.^7).*(t-25e-9).^3.*exp(-(t-25e-9).^2./pw.^2)+ ...
(16.*(t-25e-9).^5./pw.^9.*exp(-(t-25e-9).^2/pw.^2)))+ ...
A*((60./pw.^5).*(t-30e-9).*exp(-(t-30e-9).^2./pw.^2)-(80./pw.^7).*(t-30e-9).^3.*exp(-(t-30e-9).^2./pw.^2)+ ...
(16.*(t-30e-9).^5./pw.^9.*exp(-(t-30e-9).^2/pw.^2)))+ ...
-A*((60./pw.^5).*(t-35e-9).*exp(-(t-35e-9).^2./pw.^2)-(80./pw.^7).*(t-35e-9).^3.*exp(-(t-35e-9).^2./pw.^2)+ ...
(16.*(t-35e-9).^5./pw.^9.*exp(-(t-35e-9).^2/pw.^2)));


%-A inverts waveform

%DPSK delayed(5e-9) fifth(ydelayed)(notice inversion from ydiff-keeps ADC
%(comparator) output correct). This is acceptable under differential
%encoding PSK phase ambiguity rules.
ydelayed=A*yd+ ...
-A*((60./pw.^5).*(t-10e-9).*exp(-(t-10e-9).^2./pw.^2)-(80./pw.^7).*(t-10e-9).^3.*exp(-(t-10e-9).^2./pw.^2)+ ...
(16.*(t-10e-9).^5./pw.^9.*exp(-(t-10e-9).^2/pw.^2)))+ ...
-A*((60./pw.^5).*(t-15e-9).*exp(-(t-15e-9).^2./pw.^2)-(80./pw.^7).*(t-15e-9).^3.*exp(-(t-15e-9).^2./pw.^2)+ ...
(16.*(t-15e-9).^5./pw.^9.*exp(-(t-15e-9).^2/pw.^2)))+ ...
-A*((60./pw.^5).*(t-20e-9).*exp(-(t-20e-9).^2./pw.^2)-(80./pw.^7).*(t-20e-9).^3.*exp(-(t-20e-9).^2./pw.^2)+ ...
(16.*(t-20e-9).^5./pw.^9.*exp(-(t-20e-9).^2/pw.^2)))+ ...
A*((60./pw.^5).*(t-25e-9).*exp(-(t-25e-9).^2./pw.^2)-(80./pw.^7).*(t-25e-9).^3.*exp(-(t-25e-9).^2./pw.^2)+ ...
(16.*(t-25e-9).^5./pw.^9.*exp(-(t-25e-9).^2/pw.^2)))+ ...
-A*((60./pw.^5).*(t-30e-9).*exp(-(t-30e-9).^2./pw.^2)-(80./pw.^7).*(t-30e-9).^3.*exp(-(t-30e-9).^2./pw.^2)+ ...
(16.*(t-30e-9).^5./pw.^9.*exp(-(t-30e-9).^2/pw.^2)))+ ...
-A*((60./pw.^5).*(t-35e-9).*exp(-(t-35e-9).^2./pw.^2)-(80./pw.^7).*(t-35e-9).^3.*exp(-(t-35e-9).^2./pw.^2)+ ...
(16.*(t-35e-9).^5./pw.^9.*exp(-(t-35e-9).^2/pw.^2)))+ ...
A*((60./pw.^5).*(t-40e-9).*exp(-(t-40e-9).^2./pw.^2)-(80./pw.^7).*(t-40e-9).^3.*exp(-(t-40e-9).^2./pw.^2)+ ...
(16.*(t-40e-9).^5./pw.^9.*exp(-(t-40e-9).^2/pw.^2)));


ydiff=ydiff+noise;%DPSK modulated fifth with noise
ydelayed=ydelayed+noise;%delayed DPSK fifth with noise


yc=ydelayed.*ydiff;
%yc=ydiff.*ydelayed;%This is where the correlation occurs in the receiver and would be the
%mixer in the receiver. You are multipling the DPSK fifth with a .5e-9 delayed replecia 
%==================================================
% FFT
%==================================================
%new FFT for DPSK modulated fifth(ydiff)
NFFYDIFF=2.^(ceil(log(length(ydiff))/log(2)));
FFTYDIFF=fft(ydiff,NFFYDIFF);%pad with zeros
NumUniquePts=ceil((NFFYDIFF+1)/2); 
FFTYDIFF=FFTYDIFF(1:NumUniquePts);
MYDIFF=abs(FFTYDIFF);
MYDIFF=MYDIFF*2;
MYDIFF(1)=MYDIFF(1)/2;
MYDIFF(length(MYDIFF))=MYDIFF(length(MYDIFF))/2;
MYDIFF=MYDIFF/length(ydiff);
f=(0:NumUniquePts-1)*2*Fn/NFFYDIFF;


%new FFT for correlated pulses(yc)
%yc is the time domain signal output of the multiplier
%(modulated times delayed) in the correlation receiver. Plots 
%in the time domain show that a simple comparator and clocked flip flop instead of a multibit high speed flash ADC 
%may be used to recover the 1001101 signal depending on integrator design and level of
%peak voltage into mixer.
NFFYC=2.^(ceil(log(length(yc))/log(2)));
FFTYC=fft(yc,NFFYC);%pad with zeros
NumUniquePts=ceil((NFFYC+1)/2); 
FFTYC=FFTYC(1:NumUniquePts);
MYC=abs(FFTYC);
MYC=MYC*2;
MYC(1)=MYC(1)/2;
MYC(length(MYC))=MYC(length(MYC))/2;
MYC=MYC/length(yc);
f=(0:NumUniquePts-1)*2*Fn/NFFYC;

%===================================================
% PLOTS
%===================================================
%plots for DPSK modulated fifth(ydiff)
figure(1)
subplot(2,2,1); plot(t,ydiff);xlabel('TIME');ylabel('AMPLITUDE');
title('DPSK modulated pulse train(01110110)');
grid on;
%axis([-1e-9,27e-9 -1 2]);
subplot(2,2,2); plot(f,MYDIFF);xlabel('FREQUENCY');ylabel('AMPLITUDE');
grid on;
%axis([0 20e9 0 .0001]);%zoom in/out
subplot(2,2,3); plot(t,ydelayed);xlabel('TIME');ylabel('AMPLITUDE');
title('DPSK delayed pulse train(01110110)');
grid on;
%axis([-1e-9,27e-9 -1 2]);