%The noise is AWGN, to use ADSLnoise, use the function ADSLnext2 instead
%of AWGN_PSD
%freq.mat, resp.mat are inputs (WE USED TA1channel) 
load resp.mat
load freq.mat

%we want row vectors
resp=transpose(resp);
freq=transpose(freq);

%Our first plot is the two sided response of the line

[resp2side,freq2side]=Show_2sided(freq,resp);
figure
subplot(2,1,1)
plot(freq2side,20*log10(abs(resp2side)))
title('magnitude of channel response (VDSL line with 2 bridge taps)')
xlabel('freq (Hz)')
ylabel('dB')
subplot(2,1,2)
plot(freq2side,angle(resp2side))
title('phase response of the channel response')
xlabel('freq (Hz)')
ylabel('theta (radians)')
%Used to check normalizations
max(20*log10(abs(resp2side)))



%to perform bandpass filtering, we must define the parameters for the 
%bandpass filter (THESE ARE USER PARAMETERS)

fc=freq(ceil(length(freq)/2));
BW=1280*(freq(100)-freq(99));
alpha=.5;
SNR=30;
sigs=-143;% 2 sided PSD in dBm/Hz

%This is filter
Ts=1/(2*BW);
[RCresp,des_freqRC]=Freq_RC(freq,fc,Ts,alpha); %generate RC pulse
RCresp=RCresp/max(abs(RCresp));	%FILTER SHOULD HAVE UNIT MAGNITUDE
[RC2side,freqRC2side]=Show_2sided(des_freqRC,RCresp);


%plot the bandpass filter
figure
subplot(2,1,1)
plot(freqRC2side,20*log10(abs(RC2side)))
title('mag of BPF (RC freq resp)')
xlabel('frequency (Hz)')
ylabel('magnitude(dB)')
V=axis;
axis([-max(freq),max(freq),V(3),V(4)])
subplot(2,1,2)
plot(freqRC2side, angle(RC2side))
V=axis;
title('phase of BPF (RC resp)')
xlabel('frequency (Hz)')
ylabel('phase(radians)')
axis([-max(freq),max(freq),V(3),V(4)])
max(20*log10(abs(RC2side)))


%Find and plot the analytic equivalents of the filter and the channel
%The convolution of these two is the "effective analytic channel"
[RCanal,freqRCanal]=Anal_Eq(freqRC2side,RC2side);
[respanal,freqanal]=Anal_Eq(freq2side,resp2side);
figure
clf
subplot(2,2,1)
plot(freqRCanal,20*log10(abs(RCanal)))
title('mag of analytic frequency response of BPF')
xlabel('frequency(Hz)')
ylabel('magnitude(dB)')
V=axis;
axis([-max(freq),max(freq),V(3),V(4)])
subplot(2,2,3)
fmax=max(freqRCanal)
fmin=min(freqRCanal)
index=(freqanal>=fmin & freqanal<=fmax);
maxanal=max(20*log10(abs(respanal(index))));
plot(freqanal(index),20*log10(abs(respanal(index))))
title('mag of analytic chan resp in relevant regions')
xlabel('frequency(Hz)')
ylabel('magnitude(dB)')
V=axis;
axis([-max(freq),max(freq),V(3),V(4)])
subplot(2,2,2)
plot(freqRCanal, angle(RCanal))
title('phase of analytic BPF')
xlabel('frequency(Hz)')
ylabel('phase(angle)')
V=axis;
axis([-max(freq),max(freq),V(3),V(4)])
subplot(2,2,4)
plot(freqanal(index), angle(respanal(index)))
title('phase of the analytic chan resp')
xlabel('frequency(Hz)')
ylabel('phase(angle)')
V=axis;
axis([-max(freq),max(freq),V(3),V(4)])


%Show the analytic equivalents in the relevant regions

figure
plot(freqRCanal,20*log10(abs(RCanal)))
hold on
plot(freqanal(index),20*log10(abs(respanal(index))))
title('plot of the analytic channel and filter response in relevant regions')
ylabel('magnitude (dB)')
xlabel('frequency (Hz)')


%Perform filtering of the channel. THIS IS ANALYTIC EQUIVALENT (6dB
%higher than positive frequency effective channel)

[Yanal,Yanalfreq]=Chan_anal(RCanal,respanal,freqRCanal,freqanal);
figure
subplot(2,1,1)
plot(Yanalfreq,20*log10(abs((Yanal))))
title('magnitude of the analytic representation of the filtered (effective) response')
xlabel('frequency(Hz)')
ylabel('magnitude(dB)')
V=axis;
axis([-max(freq),max(freq),V(3),V(4)])
subplot(2,1,2)
plot(Yanalfreq,angle(Yanal))
title('phase of the analytic representaiton of the effective response')
xlabel('frequency(Hz)')
ylabel('phase(radian)')
V=axis;
axis([-max(freq),max(freq),V(3),V(4)])
max(20*log10(abs(Yanal)))



%Find the baseband equivalent (H_A(w+wc)) and divide by 2 since 
%Ybb(w)=.5*Hbb(w)*Xbb(w), so we can think of .5Hbb(w) filtering Xbb(w)
%Make a plot

[Ybb,Ybbfreq]=Chan_bbout(Yanalfreq,Yanal);
Ybb=Ybb/2;
figure
subplot(2,1,1)
plot(Ybbfreq,20*log10(abs(Ybb)))
title('magnitude of 1/2 baseband representation of the effective response')
xlabel('frequency (Hz)')
ylabel('magnitude (dB)')
V=axis;
axis([-max(freq),max(freq),V(3),V(4)])
subplot(2,1,2)
plot(Ybbfreq,angle(Ybb))
title('phase of 1/2 baseband representation of effective response')
xlabel('frequency(Hz)')
ylabel('phase(radians)')
V=axis;
axis([-max(freq),max(freq),V(3),V(4)])
max(20*log10(abs(Ybb)))

%Find time response and plot
[time,timeresp]=Time_Resp(Ybbfreq,Ybb);
figure
subplot(211)
plot(time,real(timeresp))
title('real part of 1/2 baseband effective channel response')
xlabel('time')
ylabel('magnitude')
subplot(212)
plot(time,imag(timeresp))
title('imaginary part of 1/2 baseband effective channel response')
xlabel('time')
ylabel('magnitude')

%Plot filtered channel response in passband (No scaling or anything)
[Yout,Youtfreq]=Chan_out(Ybbfreq,Ybb,fc);
figure
subplot(2,1,1)
plot(Youtfreq,20*log10(abs(Yout)))
title('magnitude of the passband representation of the effective response')
xlabel('frequency')
ylabel('magnitude')
V=axis;
axis([-max(freq),max(freq),V(3),V(4)])
subplot(2,1,2)
plot(Youtfreq,angle(Yout))
title('phase of the passband representation of the effective response')
xlabel('frequency')
ylabel('phase')
V=axis;
axis([-max(freq),max(freq),V(3),V(4)])






%Noise processes
%This example use AWGN

sig_squared=10^(sigs/10)%energy per dimension is 1
[NPSD,freq]=AWGN_PSD(sig_squared,freq);


%plot PSD
figure
[resp2PSD,freq2PSD]=Show_2sided(freq,NPSD);
plot(freq,10*log10(2*NPSD))	%one-sided PSD
hold on
[NPSD2side,freqqer]=Show_2sided(freq,NPSD);
plot(freqqer,10*log10(NPSD2side))
V=axis;
axis([V(1) V(2) V(3)-5 V(4)+5])
title('one-sided PSD and two sided PSD')
xlabel('frequency (Hz)')
ylabel('dBm')

%DO SCALING:  multiplying by sqrt(1/2) = 
%	      change in PSD of 1/2
NPSD=Cioffi_constant(Cioffi_constant(NPSD)); 



%Filter the Noise 
[Nout,freqout]=Noise_Filt(NPSD,RCresp,freq,des_freqRC);

%Assuming above was 2-sided input, we show the two sided output
[resp2PSD,freq2PSD]=Show_2sided(freqout,Nout);
figure
clf
plot(freq2PSD,10*log10(abs(resp2PSD)))
hold on
plot(freqout,10*log10(2*abs(Nout)))
title('magnitude of the power spectrum of the scaled noise through filter')
xlabel('frequency')
ylabel('dBm')
V=axis;
axis([-max(freq),max(freq),V(3),V(4)])


%Show PSD of analytic equivalent of noise (see B.14)

[respPSDanal,freqPSDanal]=Anal_Eq(freq2PSD,resp2PSD);
respPSDanal=2*respPSDanal;	%eqn (B.14)
figure
clf
plot(freqPSDanal,10*log10(abs(respPSDanal)))
title('magnitude of the power spectrum of the analytic filtered noise signal noise')
xlabel('frequency')
ylabel('dBm')
V=axis;
axis([-max(freq),max(freq),V(3),V(4)])

%find noise in relevant region (this is a hack to use earlier programs)

dummy=2*ones(1,(length(freqRCanal)));
[PSDanal,PSDanalfreq]=Chan_anal(dummy,respPSDanal,freqRCanal,freqPSDanal);
figure
clf
plot(freqRC2side,20*log10(abs((RC2side))))
title('magnitude of filter response')
xlabel('frequency')
ylabel('dB')
V=axis;
axis([-max(freq),max(freq),V(3),V(4)])


%Show Baseband equivalent noise spectrum
[PSDbb,PSDbbfreq]=Chan_bbout(PSDanalfreq,PSDanal);
figure
clf
plot(PSDbbfreq,10*log10(abs(PSDbb)))
title('magnitude of the power spectrum of the baseband version of the filtered noise')
xlabel('frequency')
ylabel('dBm')
V=axis;
axis([-max(freq),max(freq),V(3),V(4)])

%Show the Output PSD which should be the same as above with Show2_sided
[PSDout,PSDfreqout]=Chan_out(PSDbbfreq,PSDbb,fc);
PSDout=PSDout/2;	%(B.25)
figure
clf
plot(PSDfreqout,10*log10(abs(PSDout)),'.')
V=axis;
title('magnitude of the passband power spectrum of the filtered noise in relevant regions')
xlabel('frequency')
ylabel('dBm')
axis([-max(freq),max(freq),V(3),V(4)])


%Find equivalent noise channel
[Equiv,Equivfreq]=EquivNoise(PSDbb,Ybb,PSDbbfreq,Ybbfreq);
figure
plot(Equivfreq,20*log10(abs(Equiv)))
title('equivalent white noise channel')
xlabel('frequency(Hz)')
ylabel('magnitude (dB)')


