# Simulation script for Controlled Envelope SSB (CESSB)
# For QEX
# by David L. Hershberger W9GR
# w9gr@arrl.net
#
tic;
# Read in the peak-limited speech that was created with SOX:
[speech,fs]=wavread('SSBaudioprocessed.wav');
# The sampling rate (48 kHz) is read from the wav file.
nspeech=length(speech);
audiopeak1=max(abs(speech));
# Create a shaped 1 kHz 1 second sinewave audio tone burst and place it at the beginning of the speech
# This tone is at the peak reference level
ftone=1000;
ttone=1;
ntsamp=round(ttone*fs);
tone=sqrt(0.5)*sin(2*pi*linspace(0,(ntsamp-1)/fs,ntsamp)*ftone);
# Now shape the tone edges with a raised cosine shape
tedge=0.05;
ntedge=round(tedge*fs);
edge=0.5-0.5*cos(pi*linspace(0.5,ntedge-0.5,ntedge)/ntedge);
tone(1:ntedge)=tone(1:ntedge).*edge;
tone(ntsamp-ntedge+1:ntsamp)=tone(ntsamp-ntedge+1:ntsamp).*fliplr(edge);
# Create silence
tsilence=0.2;
ntsilence=round(tsilence*fs);
# Put silence and tone at beginning of the speech
audio=[zeros(ntsilence,1);tone.';zeros(ntsilence,1);speech]*sqrt(2);
audioclipped=[zeros(ntsilence,1);tone.';zeros(ntsilence,1);max(-sqrt(0.5),min(2*speech,sqrt(0.5)))]*sqrt(2);
ns=length(audio);

# Write the peak limited audio with the reference level tone to a wav file
# This audio file WILL create overshoot in SSB modulators
peaklimitedaudio=audio*sqrt(0.5);
wavwrite(peaklimitedaudio,fs,'peak-limited-audio.wav');


# The next block of code converts the peak limited speech (plus the reference tone)
# into CESSB baseband signals. It begins with a Weaver SSB modulator.

# Make Weaver LPF to create SSB
# The number of coefficients is large because
# the sampling rate is high - 48 kHz.
# For the same shape factor, the number of coefficients
# will go down in proportion to the reduced sampling rate.
nweaver=1820;
fnyq=fs/2;
fcut=1350;
fstop=1450;
fremez=[0,fcut/fnyq,fstop/fnyq,1];
aremez=[1,1,0,0];
weaverlpf=remez(nweaver,fremez,aremez);
nweaver=length(weaverlpf);
# Sweep and plot the resulting filter
nw=4096;
f1=0;
f2=fnyq*nw/(nw+1);
w1=2*pi*f1/fs;
w2=2*pi*f2/fs;
fweaver=linspace(f1,f2,nw);
wweaver=2*pi*fweaver/fs;
wweaverresp=freqz(weaverlpf,1,wweaver);
figure(1);
plot(fweaver,20.0*log10(abs(wweaverresp)));
title("Weaver Filter Response");
xlabel("Frequency (Hertz)");
ylabel("Response (decibels)");
axis([min(fweaver) max(fweaver) -80 5]);
grid on;
saveas(1,"Weaver-Filter-Response.png");
# Plot passband
[~,n2]=min(abs(fweaver-fcut));
figure(2);
plot(fweaver(1:n2),20*log10(abs(wweaverresp(1:n2))));
title("Weaver Filter Passband");
xlabel("Frequency (Hertz)");
ylabel("Response (decibels)");
grid on;
saveas(2,"Weaver-Filter-Passband.png");
# Set the Weaver folding frequency to be the arithmetic mean
# of the low and high frequency cutoffs
ffold=0.5*(300+3000);
time=linspace(0,(ns-1)/fs,ns).';
# Create the Weaver folding carrier
wcarrier=exp(i*time*2*pi*ffold);
# Now do the frequency folding and the Weaver filtering in the next statement
weaveriq=filter(weaverlpf,1,audio.*wcarrier);
# Compute the envelope function
envssb=2*abs(weaveriq);
# Show the envelope overshoot of the SSB modulator
figure(3);
plot(time,envssb,"linewidth",2);
axis([0,ceil(time(ns)),0,ceil(max(envssb)*20)/20]);
xlabel("Time (seconds)");
ylabel("Normalized Envelope Amplitude");
grid on;
titl=sprintf("SSB Envelope from Peak Limited Audio %0.2f%% Overshoot",100*(max(envssb)-1));
title(titl);
saveas(3,"SSB-envelope.png")

# Now do baseband envelope clipping
# Create the divider signal for baseband envelope clipping
bbenvclipdenom=max(1,envssb);
# Form the unfiltered, baseband envelope clipped signal
bbenvclip=weaveriq./bbenvclipdenom;
# bbenvclipfil is the output of the conventional "baseband envelope clipper"
# Now do the filtering to remove out of band distortion products
bbenvclipfil=filter(weaverlpf,1,bbenvclip);
# Calculate the envelope of the baseband clipped and filtered SSB signal
envbbenvclip=2*abs(bbenvclipfil);
# Now show the envelope function of the baseband envelope clipper
figure(4)
plot(time,envbbenvclip,"linewidth",2);
axis([0,ceil(time(ns)),0,ceil(max(envssb)*20)/20]);
xlabel("Time (seconds)");
ylabel("Normalized Envelope Amplitude");
grid on;
titl=sprintf("SSB Envelope After Baseband Clipping %0.2f%% Overshoot",100*(max(envbbenvclip)-1));
title(titl);
saveas(4,"SSB-clipped-envelope.png");

# Now take the output of the baseband envelope clipper and
# do the "more than clipping" to make it into CESSB:

# Make second slightly wider Weaver LPF
nweaver2=1730;
fcut2=1450;
fstop2=1550;
fremez2=[0,fcut2/fnyq,fstop2/fnyq,1];
weaverlpf2=remez(nweaver2,fremez2,aremez);
nweaver2=length(weaverlpf2);
# Sweep and plot the resulting filter
wweaverresp2=freqz(weaverlpf2,1,wweaver);
figure(5);
plot(fweaver,20.0*log10(abs(wweaverresp2)));
title("Second Weaver Filter Response");
xlabel("Frequency (Hertz)");
ylabel("Response (decibels)");
axis([min(fweaver) max(fweaver) -80 5]);
grid on;
saveas(5,"Second-Weaver-Filter-Response.png");

# env=2*abs(bbenvclipfil);
# overshootcomp is the nonlinearly overshoot-compensated signal before final filtering
lenv=length(envbbenvclip);
# Do the peak stretching
lhold=2;
envhold=envbbenvclip;
for ihold=-lhold:1:-1
  envhold=max(envhold,[zeros(-ihold,1);envbbenvclip(1:lenv+ihold)]);
endfor
for ihold=1:1:lhold
  envhold=max(envhold,[envbbenvclip(1+ihold:lenv); zeros(ihold,1)]);
endfor
# Form the clippings signal
clippings=max(1,envhold)-1;
# Set the gain factor
x1=2.0;
# Create the divider signal
overshootdiv=1+x1*clippings;
# Do the overshoot compensation
overshootcomp=bbenvclipfil./overshootdiv;
# Do the final filtering to produce the system output
cessb=filter(weaverlpf2,1,overshootcomp);
# Calculate the CESSB envelope function
envcessb=2*abs(cessb);
#

# Plot the CESSB envelope
figure(6);
plot(time,envcessb,"linewidth",2);
axis([0,ceil(time(ns)),0,ceil(max(envssb)*20)/20]);
xlabel ("Time (seconds)");
ylabel ("CESSB Envelope");
titl=sprintf("CESSB Envelope %0.2f%% Overshoot",100*(max(envcessb)-1));
title(titl);
grid on;
saveas(6,"CESSB-Envelope.png");

# Now that we have the complex cessb signal, make a real signal out of it
# Choose a random phase in degrees
cessbphase=37.5;
# The complex cessb vector is in the Weaver domain, so shift
# it back to normal frequencies while converting it to real audio
cessbaudio=2*real(cessb.*conj(wcarrier)*exp(1i*cessbphase*pi/180.));
# Now write the CESSB audio to a WAV file, so that it can be used
# to test radios for "CESSB Ready" status
# First scale the audio to -3 dB from full scale
afpeak=max(abs(cessbaudio));
cessbreadytestaudio=sqrt(0.5)*cessbaudio/afpeak;
wavwrite(real(cessbreadytestaudio),fs,'CESSB-ready-test-audio.wav');


# Make a linear phase SSB filter to produce filter method USB at 12 kHz
nssb=1022;
fedgesssb=[0 12000 12150 15150 15300 24000]*2/fs;
aedges=[0 0 1 1 0 0];
wtssb=[32 1 32];
ssbfilter=remez(nssb,fedgesssb,aedges,wtssb);
nsweepssb=5001;
fsweepssb=linspace(11000,16000,nsweepssb);
hssb=freqz(ssbfilter,1,fsweepssb*2*pi/fs);
figure(7);
plot(fsweepssb,20*log10(abs(hssb)));
title("Linear Phase SSB BPF");
xlabel("Frequency (Hertz)");
ylabel("Response (decibels)");
axis([11000 16000 -90 1]);
grid on;
saveas(7,"Linear-Phase-SSB-Response.png");

# Now set the carrier frequency
fc=12000;
nsamp=length(cessbaudio);
t=(linspace(0,nsamp-1,nsamp)/fs).';
# Do the filter method modulation and filtering in the next statement
filterssb=2*filter(ssbfilter,1,sin(2*pi*fc*t).*cessbaudio);
# Now produce a time domain plot
figure(8)
plot(t,filterssb);
# This plot is used in the article text
title("Linear Phase BPF SSB Time Domain");
xlabel("Time (seconds)");
ylabel("Amplitude");
grid on;
saveas(8,"Linear-Phase-BPF-SSB.png");

# Create the spectrum of the filter method SSB signal
nfft=4096;
[specfilterssb,fspecfilterssb]=pwelch(filterssb,[],[],nfft,fs,"half");
figure(9);
plot(fspecfilterssb,10*log10(specfilterssb));
title("Linear Phase BPF SSB Spectrum");
xlabel("Frequency (Hertz)");
ylabel("Magnitude (decibels)");
grid on;
saveas(9,"Linear-Phase-BPF-SSB-Spectrum.png");

# Now do a linear phase Hilbert transform modulator
# Load the Hilbert transform coefficients
load("hilbertssb.txt");
nfull=2400;
ffull=0.5*fs*linspace(0.5,nfull-0.5,nfull)/nfull;
hhilbert=freqz(hilbertssb,1,ffull*2*pi/fs);

# Do the Hilbert transform filtering
# The zero padding and truncation is to produce "zero delay" filtering
nzeros=floor((length(hilbertssb)-1)/2);
ssbhilbert=filter(hilbertssb,1,[cessbaudio.',zeros(1,nzeros)]).';
ssbhilbert=ssbhilbert(nzeros+1:length(cessbaudio)+nzeros);
# Now form the complex baseband signal
ssbhilbert=cessbaudio-1i*ssbhilbert;
# Modulate it onto a complex carrier, and take the real part
ssbhilbert=real(ssbhilbert.*exp(1i*2*pi*fc*t));

# Create a time domain plot
figure(10)
plot(t,ssbhilbert);
# This plot is used in the article text
title("Hilbert Transform SSB Time Domain");
xlabel("Time (seconds)");
ylabel("Amplitude");
grid on;
saveas(10,"Hilbert-Transform-SSB.png");

# Now create a spectrum
[specssbhilbert,fspecssbhilbert]=pwelch(ssbhilbert,[],[],nfft,fs,"half");
figure(11)
plot(fspecssbhilbert,10*log10(specssbhilbert));
title("Hilbert Transform SSB Spectrum");
xlabel("Frequency (Hertz)");
ylabel("Magnitude (decibels)");
grid on;
saveas(11,"Linear-Phase-Hilbert-SSB-Spectrum.png");

# Now do a linear phase Weaver modulator
# First create a Weaver lowpass filter
fweaverbw=(3150-150)/2;
nweaver3=2047;
fedgeweaver=[0 fweaverbw fweaverbw+100 24000]*2/fs;
aweaver=[1 1 0 0];
wtweaver=[1 16];
weaverfilt=remez(nweaver3,fedgeweaver,aweaver,wtweaver);
# Sweep and plot the filter just created
nsweepweaver=3001;
fsweepweaver=10*linspace(0,nsweepweaver-1,nsweepweaver);
hweaver=freqz(weaverfilt,1,fsweepweaver*2*pi/fs);
figure(12);
plot(fsweepweaver,20*log10(abs(hweaver)));
#axis([0 2500 -.25 .25]);
axis([0 2500 -90 1]);
title("Weaver Filter Response");
xlabel("Frequency (Hertz)");
ylabel("Magnitude (decibels)");
grid on;
saveas(12,"Weaver-LPF-Response.png");

# Now do the frequency folding and filtering
ffold=(150+3150)/2;
weaveriq=filter(weaverfilt,1,exp(-1i*2*pi*ffold*t).*cessbaudio);
# plot(abs(weaveriq));
# Now modulate it onto a commplex carrier, and then take the real part
weaverssb=2*real(weaveriq.*exp(1i*2*pi*(fc+ffold)*t));
# Create the spectrum plot
[specweaverssb,fspecweaverssb]=pwelch(weaverssb,[],[],nfft,fs,"half");
figure(13);
plot(fspecweaverssb,10*log10(specweaverssb));
title("Weaver Modulator SSB Spectrum");
xlabel("Frequency (Hertz)");
ylabel("Magnitude (decibels)");
grid on;
saveas(13,"Weaver-Modulator-SSB-Spectrum.png");

# Now create the  time domain plot
figure(14);
plot(t,weaverssb);
# This plot is used in the article text
title("Weaver Modulator SSB Time Domain");
xlabel("Time (seconds)");
ylabel("Amplitude");
grid on;
saveas(14,"Weaver-Modulator-SSB.png");

# Now show what happens with inappropriate (nonlinear phase) SSB modulators

# Make a nonlinear phase elliptic SSB bandpass filter
# This will simulate what a conventional crystal or mechanical filter will do
efedges=2*[12150 15150]/fs;
Rp=0.01;
Rs=80;
[eb,ea]=ellip(13,Rp,Rs,efedges);
# Now sweep the filter just created
nsweepssb=501;
fsweepssb=linspace(11000,16000,nsweepssb);
hessb=freqz(eb,ea,fsweepssb*2*pi/fs);
# Approximate group delay
gdessb=-diff(unwrap(angle(hessb)))./(2*pi*diff(fsweepssb));
f1=12100;
f2=15200;
[~,n1]=min(abs(fsweepssb-f1));
[~,n2]=min(abs(fsweepssb-f2));
# The plotyy figure will be generated at the end of this script (figure 15)
# It will show dB amplitude response and group delay of the elliptic SSB filter

# Now modulate and filter to create SSB
efilterssb=2*filter(eb,ea,sin(2*pi*fc*t).*cessbaudio);
# Find the envelope peak
efilterpeak=max(abs(efilterssb));
fprintf("Nonlinear Phase Elliptic Filter SSB Envelope Peak %0.2f%% Overshoot\n",100*(efilterpeak-1));
# Create a time domain plot
figure(16);
plot(t,efilterssb);
# This plot is used in the article text
titl=sprintf("Nonlinear Phase Elliptic Filter SSB Envelope Peak %0.2f%% Overshoot\n",100*(efilterpeak-1));
title(titl);
xlabel("Time (seconds)");
ylabel("Amplitude");
grid on;
saveas(16,"Elliptic-Filter-Modulator-SSB.png");

# Create a spectral plot
nfft=4096;
[specefilterssb,fspecefilterssb]=pwelch(efilterssb,[],[],nfft,fs,"half");
figure(17);
plot(fspecefilterssb,10*log10(specefilterssb));
title("Elliptic Filter Modulator SSB Spectrum");
xlabel("Frequency (Hertz)");
ylabel("Magnitude (decibels)");
grid on;
saveas(17,"Elliptic-Filter-SSB-Spectrum.png");

# Do a nonlinear-phase phasing type modulator
# Load in the set II coefficients from Theodor A. Prosch, DL8PT
# "A Minimalist Approximation of the Hilbert Transform," QEX,
# September-October 2012, page 25
# load "phasedifference.txt";
a2=[ -0.25038357; -0.640879872; -0.863959342; -0.955490879; -0.993211654];
b2=[ -0.07191473; -0.459269214; -0.774608418; -0.920486871; -0.977660080];
# Convolve the biquad coefficients together to make a complete IIR transfer function
# for each of I and Q
a=conv([a2(1),0,1],[a2(2);0;1]);
a=conv(a,[a2(3);0;1]);
a=conv(a,[a2(4);0;1]);
a=conv(a,[a2(5);0;1]);
# Add in unit delay to the a transfer function for the a coefficients (see Prosch article)
a=conv(a,[0 1]);
# Now convolve the b coefficients
b=conv([b2(1),0,1],[b2(2);0;1]);
b=conv(b,[b2(3);0;1]);
b=conv(b,[b2(4);0;1]);
b=conv(b,[b2(5);0;1]);
# Sweep the a and b networks
nsweep=480;
fsweep=0.5*fs*linspace(0.5,nsweep-0.5,nsweep)/nsweep;
ha=freqz(a,flipud(a),fsweep,fs);
hb=freqz(b,flipud(b),fsweep,fs);
# Calculate group delay of a network
gdphdiff=-diff(unwrap(angle(ha)))./(2*pi*diff(fsweep));
# The plotyy figure will be generated at the end of this script (figure 18)
# It will show deviation from 90 degrees phase difference, and overall group delay
#
raddeg=180/pi;
# Do the filtering of the input audio and create complex audio output
cpdphasingaudio=filter(a,flipud(a),cessbaudio)-1i*filter(b,flipud(b),cessbaudio);
# Find the peak envelope
pdphasingpeak=max(abs(cpdphasingaudio));
fprintf("Nonlinear Phase Phasing Method SSB Envelope Peak %0.2f%% Overshoot\n",100*(pdphasingpeak-1));
# Modulate onto a carrier
pdphasingssb=real(exp(1i*2*pi*fc*t).*cpdphasingaudio);
# Produce a time domain plot
figure(19);
plot(t,pdphasingssb);
# This plot is used in the article text
titl=sprintf("Nonlinear Phase Phasing Method SSB Envelope Peak %0.2f%% Overshoot\n",100*(pdphasingpeak-1));
title(titl);
xlabel("Time (seconds)");
ylabel("Amplitude");
grid on;
saveas(19,"Phase-Difference-Network-Modulator-SSB.png");
# Create a spectral plot
[specpdphasingssb,fspecpdphasingssb]=pwelch(pdphasingssb,[],[],nfft,fs,"half");
figure(20);
plot(fspecpdphasingssb,10*log10(specpdphasingssb));
title("Phase Difference Network Modulator SSB Spectrum");
xlabel("Frequency (Hertz)");
ylabel("Magnitude (decibels)");
grid on;
saveas(20,"Phase-Difference-Network-Modulator-Spectrum.png");

# Now put the processed (but NOT CESSB) audio into the elliptic BPF SSB modulator
# Now modulate and filter to create SSB
efilterssb2=2*filter(eb,ea,sin(2*pi*fc*t).*audio);
# Find the envelope peak
efilterpeak2=max(abs(efilterssb2));
fprintf("Nonlinear Phase Elliptic Filter SSB with non-CESSB audio Envelope Peak %0.2f%% Overshoot\n",100*(efilterpeak2-1));
# Create a time domain plot
figure(21);
plot(t,efilterssb2);
# This plot is used in the article text
titl=sprintf("Nonlinear Phase Elliptic Filter SSB with non-CESSB audio Envelope Peak %0.2f%% Overshoot\n",100*(efilterpeak2-1));
title(titl);
xlabel("Time (seconds)");
ylabel("Amplitude");
grid on;
saveas(21,"Elliptic-Filter-Modulator-SSB-2.png");
figure(22)
plot(t,efilterssb2,'r',t,efilterssb,'b');
# This plot is used in the article text
title("SSB Elliptic Filter Modulator with Peak Limited Audio and CESSB Audio Inputs");
xlabel("Time (seconds)");
ylabel("Amplitude");
grid on;
saveas(22,"Elliptic-Filter-Modulators-SSB.png");
fprintf("Nonlinear Phase Filter SSB Envelope Peaks %0.2f%% Overshoot %0.2f%% Overshoot\n",100*(efilterpeak-1),100*(efilterpeak2-1));


# Do the double y axis plots (plotyy):

figure(15);
ax=plotyy(fsweepssb,20*log10(abs(hessb)),fsweepssb(n1:n2),gdessb(n1:n2));
title("Amplitude and Group Delay Response of Elliptic BPF");
xlabel("Frequency (Hertz)");
ylabel(ax(1),"Amplitude (decibels)");
ylabel(ax(2),"Group Delay (seconds)");
axis(ax(1),[11000 16000 -100 5]);
grid on;
saveas(15,"Elliptic-BPF-Response.png");

figure(18)
ax=plotyy(fsweep,90+raddeg*(unwrap(angle(ha))-unwrap(angle(hb))),fsweep(1:nsweep-1),gdphdiff);
axis(ax(1),[0 fs/2 -0.2 0.2]);
title("Phase and Group Delay Response of Phase Difference Networks");
xlabel("Frequency (Hertz)");
ylabel(ax(1),"Phase Error (degrees)");
ylabel(ax(2),"Group Delay (seconds)");
saveas(18,"Phase-Difference-Network-Response.png");


toc;


