%Weaver method SSB generation
function WeaverSSBgen()

%% message

%Fs=4*44.1e3; %give some headroom for higher carrier frequencies
%BW=22.05e3;
Fs=50e3; %[Hz]
BW=4500; %[Hz]

MaxFreq=Fs/2;

%% voice mod
%t=0:1/Fs:9.355895;
%x=upsample(wavread('voice.wav').',4);

%fm=100; %[Hz]
fc=10000; %[Hz]

%% tone mod
%t=0:1/Fs:.025;
%x(1,:)=cos(2*pi*fm*t)+cos(4*pi*fm*t); %message m(t)

%sinc^2 mod
t=-.1:1/Fs:.1;
B=2e3; %[Hz]
x(1,:)=pi*B*sinc(2*pi*B*t).^2; %message m(t)

%% Filter setup

nFIR = 1024; % Pick a length long enough to give desired transition band.
WnFIRt = (BW/2)/(Fs/2); % NORMALIZED FREQ %set to be 1/2 BW for Weaver
WnFIRr = BW/(Fs/2); % NORMALIZED FREQ
TDSamplesT = fir1(nFIR-1,WnFIRt); %Hamming window
TDSamplesR = fir1(nFIR-1,WnFIRr); %Hamming window

N=2.^nextpow2(length(x)+nFIR); %used for all FFTs

figure
plot(t,x)
title('Time-domain message signal')
xlabel('Time [sec]')
ylabel('Amplitude [dimensionless]') %diag
%% Loop for USB/LSB cases


for flip=0:1

SSB=flip;  %Enter 1 for LSB and 0 for USB;

%% message freq. shift oscillators

LOhalf0=2*cos(2*pi*BW/2*t-0); %half of the maximum BW
LOhalf90=2*cos(2*pi*BW/2*t-pi/2); %half of the maximum BW

%% SSB TX L.O.

LOtx0(1,:)=2*cos(2*pi*(fc+(-1)^SSB*BW/2)*t-0); %Transmit L.O.
LOtx90(1,:)=2*cos(2*pi*(fc+(-1)^SSB*BW/2)*t-pi/2); %Transmit L.O.

%% RX L.O.
LOrx(1,:)=2*cos(2*pi*(fc)*t+0);
%LOrx(2,:)=2*cos(pi*30e3*t+0); % add additional lines for off-freq receive



%% LPF for input file (to be sure it is appropriate for the up-conversion)
%xfilt=x; %No baseband filtering

xfilt =  myfft(x,0,TDSamplesR,N,Fs,1);
xfilt = xfilt(nFIR/2:nFIR/2+length(x)-1);
title1=('Spectrum of Message Signal');
myLINfft(xfilt,0,N,Fs,title1,1,MaxFreq);  %diag
%% upconvert signals
xupconv0=xfilt.*LOhalf0;
xupconv90=xfilt.*LOhalf90;

title1=('upconverted signal pre-LPF');
mydBfft(xupconv0,xupconv90,N,Fs,title1,2,MaxFreq);  %diag

[firIFt(1,:), firIFt(2,:)] =  myfft(xupconv0,xupconv90,TDSamplesT,N,Fs,2); %filter TX carrier signal

firIFtClip(1,:) = firIFt(1,nFIR/2:nFIR/2+length(x)-1);
firIFtClip(2,:) = firIFt(2,nFIR/2:nFIR/2+length(x)-1);

title1=('upconverted signal post-LPF');
mydBfft(firIFtClip(1,:),firIFtClip(2,:),N,Fs,title1,2,MaxFreq);  %diag

%% upconvert signals to IF
xIF0=firIFtClip(1,1:length(x)).*LOtx0;
xIF90=firIFtClip(2,1:length(x)).*LOtx90;

firTDt=1/2*(xIF0+(-1)^SSB*xIF90);
title1=('TX signal output');
mydBfft(firTDt,0,N,Fs,title1,1,MaxFreq);
myLINfft(firTDt,0,N,Fs,title1,1,MaxFreq);

figure
plot(t,firTDt)
title('Time-domain SSB signal')
xlabel('Time [sec]')
ylabel('Amplitude [dimensionless]') %diag

TXenv(flip+1,:)=abs(hilbert(firTDt(1,:)));
figure
plot(t,TXenv(flip+1,:))
title('TX Envelope of SSB signal')
xlabel('Time [sec]')
ylabel('Amplitude [dimensionless]') %diag

%% Receive signals

for i=1:length(LOrx(:,1))

    xrec(i,:)=firTDt(1,1:length(x)).*LOrx(i,:);

end

title1=('Spectrum of Demodulated Signal pre-LPF');
%mydBfft(xrec(1,:),0,N,Fs,title1,1,MaxFreq);  %diag

firTDr =  myfft(xrec(1,:),0,TDSamplesR,N,Fs,1); %for one RX freq

%for 2 RX freq
%[firTDr(1,:) firTDr(2,:)] =  myfft(xrec(1,:),xrec(2,:),TDSamplesR,N,Fs,2);


firTDr = firTDr(1,nFIR/2:nFIR/2+length(x)-1);

title1=('Spectrum of Demodulated Signal post-LPF');
mydBfft(firTDr(1,:),0,N,Fs,title1,1,MaxFreq);  %diag


%% TD demod plots
figure
plot(t,firTDr(1,:))
title('Time-domain RX demod signal'),xlabel('Time [sec]'),ylabel('Amplitude [dimensionless]') %diag

%for 2 rx freq
%subplot(2,1,2),plot(t,firTDr(2,nFIR/2:end-nFIR/2))
%axis([0,t(end),-1,1])

%% envelope comparison
%RXenv=firTDr(1,:).*(firTDr(1,:)>0); %half-wave rectifier
RXenv(flip+1,:)=abs(hilbert(firTDr(1,:)));
figure
plot(t,RXenv(flip+1,:))
title('RX Envelope of SSB signal')
xlabel('Time [sec]')
ylabel('Amplitude [dimensionless]') %diag
end
figure
plot(t,RXenv(1,:)-RXenv(2,:))
title('Comparison of RX LSB/USB envelopes (subtraction)')
xlabel('Time [sec]')
ylabel('Amplitude [dimensionless]') %diag

figure
plot(t,TXenv(1,:)-TXenv(2,:))
title('Comparison of TX LSB/USB envelopes (subtraction)')
xlabel('Time [sec]')
ylabel('Amplitude [dimensionless]') %diag

end