Main Content

Spectrum Sharing Using Spectrum Sensing and Waveform Notching

This example illustrates the implementation and analysis of a technique for reducing mutual interference between non-cooperative radar and communication systems by inserting notches into the transmitted radar waveform. The frequency spectrum is analyzed by the radar system to determine regions occupied by a communications waveform, and a notched radar waveform that minimizes the power contained in those occupied regions is generated. This type of spectrum sharing approach enables coexistence of radar and communications systems operating in the same region of the frequency spectrum.

Introduction

Coexistence between radar and communications systems continues to be an active area of research as the frequency spectrums that these systems operate in become more crowded, particularly with the proliferation of 5G communications networks. Non-cooperative coexistence is an approach in which a radar and communication system are operating in the same frequency spectrum independently of one another and without any data exchange between systems.

The coexistence strategy explored in this example involves a radar and communications system operating in the same frequency band (a 100 MHz band around 3.5 GHz). The spectrum sharing behavior is realized by a cognitive radar system that adapts its waveform to avoid spectral overlap with the communications signal. The communications system operates freely and does not vary its transmission scheme based on the presence of the radar. Meanwhile, the radar detects the spectral location of the communications signal and inserts frequency nulls (referred to as "notches") into its transmit waveform, thereby reducing the mutual interference between the systems.

In this example, we focus on the signal level implementation of a spectrum sensing and waveform generation technique that make up the basis for this proposed coexistence strategy. Specifically, we implement the following two step process for generating a radar waveform with notches in the portion of the frequency spectrum occupied by the communications signal [1]:

  1. Determine the optimal frequency region for radar waveform transmission in the presence of an Orthogonal Frequency-Division Multiplexing (OFDM) communications waveform based on a Fast Spectrum Sensing (FSS) algorithm [5].

  2. Create and analyze a Pseudo-Random Optimized Frequency Modulation (PRO-FM) pulsed waveform with notches (areas of reduced spectrum power) in the region occupied by the communications waveform [2].

Once the process for spectrum sensing and waveform generation is examined, a simulation is run to determine the performance benefit of using a notched radar waveform in the presence of a communications signal.

Spectrum Sensing

In the coexistence scheme described in this example, a radar system minimizes the impact that its transmit waveform has on a communication system while simultaneously improving radar performance by notching the radar transmit waveform. In order for such a coexistence scheme to be effective, the radar system must continuously and quickly (on the order of a pulse repetition interval) determine an unoccupied region of the spectrum for transmission.

A popular FSS algorithm is implemented in this example [5]. This approach groups the spectrum into regions of high and low power interference plus noise, and then selects the transmit region by maximizing an objective function that balances the trade-off between the total bandwidth of the radar waveform with the anticipated Signal to Interference plus Noise Ratio (SINR) based on measurements of the spectrum that contains a communication waveform.

Generate Communications Signal

The communications signal occupying the frequency spectrum is made up of two non-overlapping OFDM transmissions. Each of these OFDM transmissions contains four carriers separated by 1 MHz (occupying a total of 4 MHz of bandwidth).

fc = 3.5e9;                            % Communications center frequency
separationFrequency = 1e6;             % 1 MHz carrier separation frequency
numCarriers = 4;                       % 4 carriers with 1 MHz separation
ofdmLocations = [fc-25e6;fc+25e6];     % One OFDM transmission at center frequency - 25 MHz, one at center frequency + 25 MHz
sampleRate = 200e6;                    % 200 MHz sample rate
numSamples = 400;                      % Collect 400 communications samples

% Helper used to generate communications signal
commsSignal = helperGenerateCommunicationsSignal(sampleRate,numCarriers,numSamples,separationFrequency,ofdmLocations,fc);

Once the communications signal has been generated, calculate and plot its power spectrum.

% Get the communications signal power spectrum
pComms = pspectrum(commsSignal,sampleRate);
pCommsDb = pow2db(pComms);

% Get the frequency samples of communications spectrum
fComms = helperSpectrumFrequencies(sampleRate,numel(pComms));

% Plot the power spectrum of the communications signal
plot(fComms/10^6,pCommsDb,'LineWidth',1);
xlabel('Frequency (MHz)'); ylabel('Power Spectrum (dB)');
title('Power Spectrum of Communications Signal');

The power spectrum of the communications signal shows two occupied frequency bands. The ideal radar waveform in the coexistence scheme described in this example would contain very little power in the regions of the frequency spectrum where the communications signals are transmitted.

Separate Spectrum into High and Low Power Meso-Bands

The first step in the FSS algorithm is to divide the spectrum into high and low power interference plus noise regions. These regions are referred to as "meso-bands". No noise is added in this example, so in this case the spectrum power is solely due to the presence of the communications signal. The algorithm for splitting the spectrum into high and low power meso-bands is:

  1. Create an alternating set of high and low power interference plus noise regions called meso-bands by comparing each sample to a user specified power threshold.

  2. For each low power region, check if it is narrower than a user specified minimum meso-band bandwidth. If a low power region has a width lower than the user specified value, merge it with the adjacent high power bands.

First, set up the power threshold and minimum bandwidth based on the values used in [1].

pCommsDbMax = max(pCommsDb);     % Get the maximum communications signal power level
pThresholdDb = pCommsDbMax-15;   % User specified power threshold is set to be 15 dB below the max communications signal level
minBandwidth = 4e6;              % Minimum meso-band bandwidth is set to be 4 MHz

A helper function is used to separate the spectrum into high and low power meso-bands using the algorithm described above.

mesoBands = helperMesoBands(pComms,pThresholdDb,minBandwidth,sampleRate);

A plot of meso-bands shows that the frequency spectrum is separated into only 5 high and low power bands. This reduces the computational complexity of determining the optimal transmission region. Grouping the frequency spectrum into high and low power meso-bands using the user specified power threshold and minimum bandwidth heuristics is what makes this spectrum sensing algorithm computationally efficient.

a = helperPlotMesoBands(pComms,mesoBands,pThresholdDb,fComms,"Merged Meso-Bands",[-50 50]);

Select Meso-Band Combination that Maximizes Objective Function

The goal of this spectrum sensing algorithm is to transmit in the frequency range that maximizes the following objective function [5]:

Objective=α*SINR+(1-α)*Bandwidth

where SINR and Bandwidth are the normalized SINR and bandwidth, and α is a user specified weighting term that balances the tradeoff between increased bandwidth and increased SINR. The following equation is used to calculate SINR for a combination of meso-bands:

SINR=τ*PsPIN

where τ is the pulse width as a fraction of the pulse repetition interval, Ps is the received signal power in a given meso-band based on the desired spectrum, and PIN is the received interference plus noise power in one pulse width in a given meso-band.

This objective function is calculated for every meso-band combination, and the combination that maximizes the objective function is selected as the radar waveform transmission region. In this case, there are only 5 meso-bands, so there are only 25-1=31 combinations of bands to calculate the objective value for.

For the purposes of this example, assume that the received radar signal power is 20 dB below the peak of the received communications signal power, the radar waveform pulse width is 5% of the pulse repetition interval (PRI), and the alpha value is set to 0.5.

peakSinrDb = -20;     % SINR of peak communications signal power vs. peak radar signal power
tau = 0.05;           % Pulse width as a fraction of pulse repitition interval
alpha = 0.5;          % Weighting factor between SINR and bandwidth

In order to determine the optimal transmission region, the desired waveform spectrum power is assumed to have a Gaussian shape. This spectral shape is desirable for radar waveforms because it results in low range sidelobes. This power spectrum is scaled to maintain the specified SINR.

% The desired waveform spectrum power has a Gaussian shape.
pWaveform = gausswin(numSamples);
pWaveformDb = pow2db(pWaveform);

% Scale the spectrum power for desired SINR
pWaveformScaledDb = pWaveformDb + (pCommsDbMax-max(pWaveformDb)) + peakSinrDb;
pWaveformScaled = db2pow(pWaveformScaledDb);

Plot the Gaussian radar waveform spectrum power along with the communications signal spectrum power used to calculate the objective function for all meso-band combinations.

% Plot the communications signal spectrum power along with the desired radar waveform
% spectrum power
helperPlotCommsWithRadarSpectrum(pComms,pWaveformScaled,'Full Radar Waveform',sampleRate,'Radar Waveform and Communications Signal Spectrum');

This plot shows that by using the full spectrum for radar transmission, there would be mutual interference between the radar waveform and the communications waveform, resulting in performance degradation for both systems.

Next, select the transmission region that maximizes the objective function. The notch frequencies are the frequencies excluded from the transmission region.

% Calculate the desired location of frequency notches in the waveform
notchFrequencies = helperTxNotchFrequencies(mesoBands,pWaveformScaled,pComms,tau,alpha,sampleRate);

Once the location of the frequency notches are determined, a desired waveform spectrum can be generated by inserting notches into the full waveform spectrum. Plot the desired notched radar waveform spectrum power along with the communications signal spectrum power.

% Generate notched spectrum
desiredRadarSpectrum = helperNotchSpectrum(pWaveformScaled,sampleRate,notchFrequencies);

% Plot communications signal spectrum along with the radar waveform spectrum desired
% for mutual interference mitigation
helperPlotCommsWithRadarSpectrum(pComms,desiredRadarSpectrum,'Desired Radar Waveform',sampleRate,'Desired Radar Waveform Spectrum based on FSS Algorithm');

This plot illustrates that using this FSS algorithm, a desired radar waveform spectrum is generated that has notches in the region of highest communications signal power. A radar waveform with this spectrum would substantially reduce mutual interference between the radar and communications systems.

Notched Waveform Design and Analysis

After the desired spectrum is determined based on the FSS algorithm, a waveform must be generated that conforms to this spectrum. In this example, we use a PRO-FM waveform [1]. This is a noise waveform that contains desirable range sidelobe characteristics and is amenable to spectral notching due to the nature of its construction. These characteristics make this a versatile waveform that is a good candidate for spectrum sharing. To generate this waveform, a seed waveform is optimized iteratively to contain a desired frequency spectrum shape and uniform amplitude [2].

Generate Initial PCFM Waveform

In the process of generating a PRO-FM waveform, you start with a seed waveform made up of data symbols with random phase. The frequency spectrum of this waveform is dictated by the symbol rate. This seed waveform is a Polyphase-Coded FM (PCFM) waveform that has random frequency values within the desired frequency band throughout the duration of the pulse [3]. This type of waveform has been shown to be amenable to spectral shaping via the PRO-FM process [2].

The transmit waveform has a design bandwidth of 100 MHz, a pulse width of 2 us, a pulse repetition interval of 80 us, and a coherent processing interval (CPI) with 200 pulses.

rng("default")            % Set random number generator for reproducibility 
bandwidth = 100e6;        % Design bandwidth is 100 MHz
pulseWidth = 2e-6;        % Pulse width is 2 us
pri = 80e-6;              % PRI is 80 us
numPulses = 200;          % 200 pulses in a CPI

% Sample at 2x our design bandwidth
sampleRate = 2*bandwidth;

The following steps describe the process of generating the PCFM waveform with random frequencies throughout the pulse duration.

  1. Randomly generate N+1 symbol phases θ,whereθ[-ππ]

  2. Calculate N phases differences αn,whereαn=θn-θn-1 for n = 1,...,N. If |αn| > π, αn=αn-2π*sgn(αn). In this equation, sgn is the sign function.

  3. Create an impulse train of phase differences using the equation n=1Nαn*δ(t-(n-1)Tp), where Tp is the time separation between impulses and δ(t) is a unit impulse at time t. This impulse train carries the random frequency information of the desired signal.

  4. Filter the impulse train to create smooth frequency transitions between samples. In this case, we use a raised-cosine filter.

  5. Integrate the filtered impulse train F(t) from 0 to t to recover the desired signal phase values. Assume θ0 as the initial condition: P(t)=0tF(t)dτ+θ0

  6. Recover the complex signal from the phase values: s(t)=eiP(t), where s(t) is the final PCFM waveform.

The following code generates the seed waveform using the algorithm described above, and plots the waveform characteristics.

% Generate the PCFM seed waveform
pcfmWaveform = helperSeedWaveform(bandwidth,pulseWidth,sampleRate,numPulses);

% Plot the waveform spectrogram for a single pulse
stft(pcfmWaveform(:,1),sampleRate); title('PCFM Seed Waveform Spectrogram');

The spectrogram of the PCFM signal shows that the frequency is randomly distributed as a function of time.

% Plot the seed waveform power spectrum averaged across pulses
pSeed = pspectrum(pcfmWaveform,sampleRate);
pSeedMean = mean(pSeed,2);
fSeed = helperSpectrumFrequencies(sampleRate,numel(pSeedMean));
plot(fSeed/10^6,pow2db(pSeedMean)); xlabel('Frequency (MHz)'); ylabel('Power Spectrum (dB)'); title('PCFM Seed Waveform Mean Spectrum Power');

The power spectrum clearly illustrates that more power is concentrated towards the center of the spectrum, although the spectral shape is not Gaussian.

% Plot the ambiguity function for a single pulse
ambgfun(pcfmWaveform(:,1),sampleRate,1/pulseWidth); title('PCFM Seed Waveform Ambiguity Function'); ylim([-3 3]); xlim([-0.3 0.3]);

Some Range-Doppler sidelobes are visible in the ambiguity function and appear to be distributed randomly, which is not unexpected due to the nature of the waveform's construction.

Generate PRO-FM Notched Waveform

The generation of the PRO-FM waveform is an iterative process where a desired frequency spectrum mask is repeatedly applied to the seed waveform signal in the frequency domain and then converted back to the time domain while maintaining constant amplitude. The following steps describe the generation of the PRO-FM waveform:

  1. Convert the previous constant amplitude time domain signal (sk) to the next time domain signal (rk+1). First convert to the frequency domain using the Fast Fourier Transform (FFT), extract phase (), apply desired spectrum (g), and perform the Inverse Fast Fourier Transform (IFFT): rk+1=ifft(g*eifft(sk))

  2. Enforce constant amplitude of time domain signal: sk+1=eirk+1

  3. Repeat steps (1) and (2) for K iterations

K is set to 500 for this example [1]. In practice, the number of iterations would be empirically determined or the iterative process could continue until the waveform characteristics meet some stopping criteria. The desired spectrum (g) can take on any shape depending on the requirements of the radar system. In this case, the notched Gaussian spectrum determined using the FSS algorithm from the previous section is used to generate the PRO-FM waveform.

iterations = 500;                           % Go through 500 waveform iterations
notchedSpectrum = desiredRadarSpectrum;     % The desired spectrum is the notched spectrum from the FSS section

The algorithm described above is used to generate the notched waveform PRO-FM waveform, and the waveform characteristics are plotted.

% Generate a PRO-FM waveform with a notched frequency spectrum
proFmNotched = helperProFm(pcfmWaveform,notchedSpectrum,iterations);

% Plot the waveform spectrogram for a single pulse
stft(proFmNotched(:,1),sampleRate); title('PRO-FM Notched Waveform Spectrogram');

The PRO-FM spectrogram still shows a randomly varying frequency throughout the duration of the pulse, although the power is more concentrated towards the center of the distribution and there are regions of low power where the frequency spectrum notches are.

% Plot the PRO-FM waveform power spectrum averaged across pulses
pProFm = pspectrum(proFmNotched,sampleRate);
pProFmMean = mean(pProFm,2);
fProFm = helperSpectrumFrequencies(sampleRate,numel(pProFmMean));
plot(fProFm/10^6,pow2db(pProFmMean)); xlabel('Frequency (MHz)'); ylabel('Power Spectrum (dB)'); title('PRO-FM Notched Waveform Mean Spectrum Power');

The power spectrum clearly illustrates the Gaussian shape of the spectrum and the presence of notches at the locations determined by FSS.

% Plot the ambiguity function for a single pulse
ambgfun(proFmNotched(:,1),sampleRate,1/pulseWidth); title('PRO-FM Notched Waveform Ambiguity Function'); ylim([-3 3]); xlim([-0.3 0.3]);

The ambiguity function shows the presence of significant Range-Doppler sidelobes. The issue of sidelobes in the waveform ambiguity function is largely mitigated by altering the transmit waveform in a CPI on a pulse to pulse basis. Using this variable waveform approach, the location of these Range-Doppler sidelobes change throughout a CPI and their impact on performance is significantly reduced [1].

Impact of Notches on Radar Performance

The previous two sections in this example showed how to determine the waveform transmission region using an FSS algorithm and how to subsequently generate a PRO-FM notched waveform using the desired transmission spectrum obtained from the FSS algorithm.

This section will analyze the impact that waveform notching has on radar performance by simulating the transmission of the notched vs. unnotched waveforms in the presence of the communications signal.

The communications signal is assumed to occupy a fixed region for the duration of a radar CPI for simplicity. This example could be adjusted so that the location of the OFDM transmission changes throughout the CPI, but for this example only the simplest case is examined. The communications signal generated in the prior section is reused here.

Generate Unnotched and Notched Waveforms

In order to analyze the impact of waveform notching on radar performance, we generate unnotched and notched PRO-FM waveforms. A new waveform is generated for each pulse in the CPI. The notched waveform has notches placed at the frequencies determined using the FSS algorithm.

% Generate notched and unnotched waveforms for a 200 pulse CPI.
unnotchedWaveformCpi = helperGenerateUnnotchedWaveform(bandwidth,pulseWidth,sampleRate,numPulses);
notchedWaveformCpi = helperGenerateNotchedWaveform(bandwidth,pulseWidth,sampleRate,numPulses,notchFrequencies);

Plot the unnotched and notched waveform spectrums alongside the communications signal spectrum to visualize the spectrum of the waveforms that are being used for transmission in this section.

% Plot communications signal as well as notched and unnotched waveform
helperPlotCommsWithRadar(commsSignal,unnotchedWaveformCpi,notchedWaveformCpi,peakSinrDb,sampleRate);

Simulate Transmission of Waveforms and Analyze Performance

In this section, we run a simulation to analyze the performance of the unnotched and notched waveforms in the absence and presence of the communications signal.

Prior to running this simulation, we create a set of transmit radar waveforms combined with communications waveforms. An alternative approach would be to add the communications signal to the entire receive pulse instead of just the transmit waveform. In this case, the communications signal is added only to the transmit waveform for simplicity to easily maintain the desired SINR on receive. Prior to combining the radar waveform with the communications waveform, this function scales the signals so that the SINR is set to the desired level.

% Add communications signal to waveforms with desired SINR
notchedWaveformCpiWithComms = helperCombineSignals(notchedWaveformCpi,commsSignal,peakSinrDb);
unnotchedWaveformCpiWithComms = helperCombineSignals(unnotchedWaveformCpi,commsSignal,peakSinrDb);

Once the radar waveforms have been generated and combined with the interfering communications waveform, set up target locations and velocities. Use 6 targets total with some moving towards the platform and some moving away from the platform to simulate a roadway.

% Set up target ranges and velocities
targetPositions = [[500;0;0],[500;0;0],[600;0;0],...
                   [600;0;0],[700;0;0],[700;0;0]];

targetVelocities = [[5;0;0],[-5;0;0],[10;0;0],...
                    [-10;0;0],[15;0;0],[-15;0;0]];

% Velocity and range limits for the Range-Doppler Plot
rLim = [450 750];
vLim = [-30 30];

After generating the waveforms and setting up the target properties, run the simulation of a single CPI and plot the Range-Doppler response of the unnotched and notched waveforms in the absence and presence of interference. For more information on simulating the transmission of custom waveforms, see Waveform Design to Improve Range Performance of an Existing System.

The first two plots show the Range-Doppler map for the notched and unnotched waveforms in the absence of the communications signal.

% Run simulation and plot Range-Doppler response for unnotched waveform
% without communications signal
[rdUnnotchedNoComms,rdRange,rdVelocity] = helperRunSimulation(unnotchedWaveformCpi,unnotchedWaveformCpi,fc,sampleRate,pri,targetPositions,targetVelocities,"Unnotched Radar Waveform, No Comms Signal",vLim,rLim);

% Run simulation and plot Range-Doppler response for notched waveform
% without communications signal
rdNotchedNoComms = helperRunSimulation(notchedWaveformCpi,notchedWaveformCpi,fc,sampleRate,pri,targetPositions,targetVelocities,"Notched Radar Waveform, No Comms Signal",vLim,rLim);

Higher range sidelobes are evident in the notched waveform when compared with the unnotched waveform. Without the presence of the communications signal, the performance of the unnotched waveform qualitatively appears to be better than the notched waveform, which is expected.

The plots below show performance of the notched and unnotched waveforms in the presence of the communications signal. This is where the performance benefits of the notched waveform are evident.

% Run simulation and plot Range-Doppler response for unnotched waveform
% with communications signal
rdUnnotchedWithComms = helperRunSimulation(unnotchedWaveformCpiWithComms,unnotchedWaveformCpi,fc,sampleRate,pri,targetPositions,targetVelocities,"Unnotched Radar Waveform With Comms Signal",vLim,rLim);

% Run simulation and plot Range-Doppler response for unnotched waveform
% with communications signal
rdNotchedWithComms = helperRunSimulation(notchedWaveformCpiWithComms,notchedWaveformCpi,fc,sampleRate,pri,targetPositions,targetVelocities,"Notched Radar Waveform With Comms Signal",vLim,rLim);

When using an unnotched waveform in the presence of the communications signal, the targets become difficult to distinguish on the Range-Doppler map. Contrast this with the notched waveform Range-Doppler map, where the targets are still clearly visible.

In order to quantify the performance difference between each test case, the average SINR of all the targets is compared. For this example, the noise level is calculated in an unoccupied region from 620 to 680 meters in range and -10 to 10 m/s in Doppler velocity.

% Use 620 to 680 m, -10 to 10 m/s Doppler, as the region to calculate noise
% power (excluding target locations)
noiseRange = [620 680];
noiseVelocity = [-10 10];

% Display a table of target SINR values as well as differences between each
% test case and the base case of unnotched waveform no communications signal.
T = helperGetSinrTable(rdUnnotchedNoComms,rdNotchedNoComms,rdUnnotchedWithComms,rdNotchedWithComms,rdRange,rdVelocity,targetPositions,targetVelocities,noiseRange,noiseVelocity);
disp(T);
                                    SINR at Target Locations (dB)    Delta with Unnotched, No Comms Signal
                                    _____________________________    _____________________________________

    Unnotched, No Comms Signal                 45.797                                     0               
    Notched, No Comms Signal                   43.724                               -2.0728               
    Unnotched, With Comms Signal               15.623                               -30.173               
    Notched, With Comms Signal                 28.518                               -17.278               

This quantitative analysis confirms what is observed qualitatively in the Range-Doppler plots - without any communications signal, the unnotched radar waveform has a marginally higher SINR than the notched radar waveform. However, when the communications signal is present, the SINR of the notched radar waveform is significantly higher than that of the unnotched radar waveform.

Using this type of notching, radar systems can coexist with communication systems while continuing to detect targets effectively.

Summary

In this example, you learned how to implement a fast spectrum sensing algorithm to determine the desirable frequency transmission region, as well a waveform notching algorithm to limit transmit frequency to the desirable region. The effectiveness of this approach was then analyzed by simulating the transmission, propagation, and reception of the custom waveform in the presence of a communications signal. This type of workflow may become more common as the need for coexistence schemes between radar and communications systems become more pressing in the future.

References

[1] Ravenscroft et al. Experimental demonstration and analysis of cognitive spectrum sensing and notching for radar. IET Radar, Sonar & Navigation. 2018.

[2] Jakabosky et al. Waveform Design and Receive Processing for Nonrecurrent Nonlinear FMCW Radar. IEEE. 2015.

[3] Blunt et al. Polyphase-Coded FM (PCFM) Radar Waveforms, Part 1: Implementation. IEEE. 2014.

[4] Higgins et al. Mitigating Interference via Spatial and Spectral Nulls. IEEE International Conference on Radar Systems. 2012.

[5] Martone et al. Spectrum Allocation for Non-Cooperative Radar Coexistence. IEEE Transactions on Aerospace and Electronic Systems. 2018.

Helper Functions

The following group of helper functions implement the spectrum sensing algorithms discussed in [5].

function mesoBands = helperMesoBands(interferenceSpectrumPower,powerThresholdDb,minBandwidth,sampleRate)
    % Get the initial meso bands in high and low power
    initialMesoBands = helperInitialMesoBands(interferenceSpectrumPower,powerThresholdDb);

    % Merge high and low power meso bands
    mesoBands = helperMergeMesoBands(initialMesoBands,minBandwidth,sampleRate);
end

function initialMesoBands = helperInitialMesoBands(spectrumPower,powerThresholdDb)
% Generate initial meso-bands from a power spectrum for a given
% threshold    

% Get the spectrum power in dB
spectrumPowerDb = pow2db(spectrumPower);
    
% Hold meso-band start/end index, and whether it is high power
initialMesoBands = [];

% Initialize loop parameters
spectrumSamples = length(spectrumPowerDb);
iStart = 1;
mesoBandIdx = 1;

for iSample = 1:spectrumSamples
    % Check if the current band start and stop powers are high
    startHighPower = spectrumPowerDb(iStart) > powerThresholdDb;
    stopHighPower = spectrumPowerDb(iSample) > powerThresholdDb;

    % If the start and stop values are in the same power band, advance the
    % stop index continue, otherwise this is the end of a meso-band.
    if startHighPower ~= stopHighPower
        initialMesoBands = helperSetMesoBands(initialMesoBands,mesoBandIdx,iStart,iSample-1,startHighPower);
        mesoBandIdx = mesoBandIdx + 1;
        iStart = iSample;
    end
end

% Set the last meso-band
initialMesoBands = helperSetMesoBands(initialMesoBands,mesoBandIdx,iStart,spectrumSamples,spectrumPowerDb(iStart) > powerThresholdDb);

end

function mergedMesoBands = helperMergeMesoBands(initialMesoBands,minBandwidth,sampleRate)
% Merge low power meso-bands for a given minimum bandwidth

% Get the number of samples in the spectrum
interferenceSamples = initialMesoBands(end).iStop-initialMesoBands(1).iStart+1;

% Calculate the minimum number of samples based on minimum bandwidth
Lmin = floor(interferenceSamples/sampleRate*minBandwidth);

% Set the initial index based on whether 1st meso-band is high power
mergedMesoBands = initialMesoBands;
if mergedMesoBands(1).HighPower
    currentIdx = 2;
else
    currentIdx = 1;
end

[~,numBands] = size(mergedMesoBands);

% Merge low power bands if they are too narrow
while currentIdx <= numBands
    currentBand = mergedMesoBands(currentIdx);
    currentBandSampleNum = currentBand.iStop - currentBand.iStart + 1;
    if currentBandSampleNum < Lmin
        [mergedMesoBands,currentIdx,numBands] = mergeSingleMesoBand(mergedMesoBands,currentIdx);
    else
        % Increment by 2 to ensure we only look at low power bands
        currentIdx = currentIdx + 2;
    end
end

% Merge the meso-band at the given index with neighbors
function [newMesoBands,currentIdx,newNumBands] = mergeSingleMesoBand(mesoBands,currentIdx)
    [~,initialNumBands] = size(mesoBands);
    
    % Determine the range of bands to combine. This is different if we
    % are on the first or the last meso-band.
    if currentIdx == 1
        newNumBands = initialNumBands-1;
        currentIdx = 2;
        combineBandRange = 1:2;
    elseif currentIdx == initialNumBands
        newNumBands = initialNumBands-1;
        combineBandRange = initialNumBands-1:initialNumBands;
    else
        newNumBands = initialNumBands-2;
        combineBandRange = currentIdx-1:currentIdx+1;
    end
    
    % Get the start and end index of the meso-bands that are being
    % combined.
    startIdx = mesoBands(combineBandRange(1)).iStart;
    endIdx = mesoBands(combineBandRange(end)).iStop;
    
    % Create the new meso-band, which is always high power
    newBand = helperCreateMesoBand(startIdx,endIdx,true);

    % Get the indexes of the prior and next bands after combining
    priorBandIdx = combineBandRange(1)-1;
    nextBandIdx = combineBandRange(end)+1;

    % Set the new meso-bands after combining the specified band
    newMesoBands = [];
    if priorBandIdx > 0
        newMesoBands = mesoBands(1:priorBandIdx);
    end

    newMesoBands = [newMesoBands,newBand];

    if nextBandIdx <= initialNumBands
        newMesoBands = [newMesoBands,mesoBands(nextBandIdx:end)];
    end  
end
end

function mesoBands = helperSetMesoBands(mesoBands,mesoBandIdx,iStart,iStop,highPower)
% Create meso-band struct row, add it to the meso-bands struct

if isempty(mesoBands)
    mesoBands = helperCreateMesoBand(iStart,iStop,highPower);
else
    mesoBands(mesoBandIdx) = helperCreateMesoBand(iStart,iStop,highPower);
end
end

function newBand = helperCreateMesoBand(startIdx,endIdx,highPower)
% Create a new meso-band

newBand = struct('iStart',startIdx,'iStop',endIdx,'HighPower',highPower);
end

function a = helperPlotMesoBands(spectrumPower,mesoBands,threshold,spectrumFrequencies,plotTitle,plotXlim)
% Plot meso-bands along with the frequency spectrum

% Offset the high and low power bands by 10 dB for display purposes
offset = 10;

% Get the spectrum power in dB
spectrumPowerDb = pow2db(spectrumPower);

% Get the start and stop idxs
threshIdxs = [mesoBands(1).iStart mesoBands(end).iStop];

% Get spectrum frequencies in MHz
spectrumFrequencies = spectrumFrequencies/1e6;

% Plot the spectrum and meso-bands
a = axes(figure);
hold(a,"on");
lThresh = plot(a,spectrumFrequencies(threshIdxs),[threshold threshold],'DisplayName','Threshold','LineStyle','--','Color','k','LineWidth',1);
lSpec = plot(a,spectrumFrequencies,spectrumPowerDb,'DisplayName','Measured Spectrum','LineWidth',1,'Color','b');
hold(a,"off");
lLow = plotMesoBandLevel(mesoBands(~[mesoBands.HighPower]),threshold,threshold-offset,spectrumFrequencies,"Low Power Meso-Band","g",a);
lHigh = plotMesoBandLevel(mesoBands([mesoBands.HighPower]),threshold,threshold+offset,spectrumFrequencies,"High Power Meso-Band","r",a);
ylabel(a,'Power (dB)')
xlabel(a, 'Frequency (MHz)')
legend(a,[lSpec,lLow,lHigh,lThresh],"Location","southeast");
title(a,plotTitle);
xlim(a,plotXlim);

% Plot high or low power meso-bands
function l = plotMesoBandLevel(mesoBands,threshold,plotPower,freq,description,color,a)
    [~,numBands] = size(mesoBands);

    hold(a,"on");
    for iBand = 1:numBands
        start = mesoBands(iBand).iStart;
        stop = mesoBands(iBand).iStop;
        range = start:stop;
        x = [start,range,stop];
        y = [threshold,repmat(plotPower,[1,length(range)]),threshold];
        l = plot(a,freq(x),y,'DisplayName',description,'Color',color,'LineWidth',2);
    end
    hold(a,"off");
end

end

function notchFrequencies = helperTxNotchFrequencies(mesoBands,signalSpectrumPower,interferenceSpectrumPower,tau,alpha,sampleRate)
% Determine notch frequency locations by maximizing the objective function
% for all combinations of messo-bands
% Determine the frequencies to notch
            
% Calculate the tx band objective values
txBandObjValues = getTxBandObjectiveValues(mesoBands,signalSpectrumPower,interferenceSpectrumPower,tau,sampleRate);

% Calculate normalized objective values for each tx band
txBandObjValuesNorm = normalizeTxBandObjectiveValues(txBandObjValues);

% Get the transmit region with the max objective function
txMesoBands = getMaxObjectiveBands(txBandObjValuesNorm,alpha);

% Get the notch frequency regions
spectrumFrequencies = helperSpectrumFrequencies(sampleRate,numel(interferenceSpectrumPower));
notchFrequencies = getNotchFrequencyFromTxRegion(mesoBands,txMesoBands,spectrumFrequencies);

function txBandObjectives = getTxBandObjectiveValues(mesoBands,signalSpectrumPower,interferenceSpectrumPower,tau,sampleRate)
    % Get all possible band combinations
    [~,numBands] = size(mesoBands);
    possibleTxBands = getAllPossibleTxBands(numBands);

    % Calculate objectives for each tx band combination
    txBandObjectives = calculateTxBandObjectiveValues(possibleTxBands,mesoBands,signalSpectrumPower,interferenceSpectrumPower,tau,sampleRate);

    function possibleBands = getAllPossibleTxBands(numBands)
        % Get every permutation of meso-band combinations
        possibleBands = struct('BandCombo',[],'SINR',[],'Bandwidth',[]);
        bandComboIdx = 1;
        
        availableBands = 1:numBands;
        
        for iBandCombos = 1:numBands
            combos = nchoosek(availableBands,iBandCombos);
            [numCombos,~] = size(combos);
            for iCombo = 1:numCombos
                possibleBands(bandComboIdx,1).BandCombo = combos(iCombo,:);
                bandComboIdx = bandComboIdx+1;
            end
        end
    end

    function txBandCalculations = calculateTxBandObjectiveValues(txBandCalculations,mesoBands,signalSpectrumPower,interferenceSpectrumPower,tau,sampleRate)
        % Calculate the objective functions for each possible meso-band combination
        
        [numPossibleBands,~] = size(txBandCalculations);

        % Get the frequency values for interference and signal spectrums
        fInterference = helperSpectrumFrequencies(sampleRate,numel(interferenceSpectrumPower));
        fStepInterference = fInterference(2)-fInterference(1);
        fSignal = helperSpectrumFrequencies(sampleRate,numel(signalSpectrumPower));
        fStepSignal = fSignal(2)-fSignal(1);
        
        % Set the sinr and bandwidth for all possible bands
        for iBandCombo = 1:numPossibleBands
            % Get the indexes of the interference and signal in this
            % meso-band
            bands = txBandCalculations(iBandCombo).BandCombo;
            [intBandIdxs,sigBandIdxs] = getBandIdxs(mesoBands,bands,fInterference,fSignal);
    
            % Get the power contained in the signal. Multiplication by tau in order to
            % capture effect of pulse width.
            mesoBandSignalPower = sum(signalSpectrumPower(sigBandIdxs))*tau*fStepSignal;
            
            % Get the power contained in the interference
            mesoBandInterferencePower = sum(interferenceSpectrumPower(intBandIdxs))*fStepInterference;
            
            % Calculate SINR and bandwidth for this mesoband
            sinr = mesoBandSignalPower / mesoBandInterferencePower;
            bandwidth = length(intBandIdxs)*fStepInterference;
    
            % Store the SINR and bandwidth
            txBandCalculations(iBandCombo).SINR = sinr;
            txBandCalculations(iBandCombo).Bandwidth = bandwidth;
        end

        function [intBandIdxs,sigBandIdxs] = getBandIdxs(mesoBands,bands,fInterference,fSignal)
        % Get band indexes from meso-bands.
            intBandIdxs = [];
            sigBandIdxs = [];
            for iBand = 1:length(bands)
                currentIntIdxs = mesoBands(bands(iBand)).iStart:mesoBands(bands(iBand)).iStop;
                currentSigIdxs = getIdxFromFreq(fInterference(currentIntIdxs(1)),fInterference(currentIntIdxs(end)),fSignal);
                intBandIdxs = [intBandIdxs,currentIntIdxs];
                sigBandIdxs = [sigBandIdxs,currentSigIdxs];
            end

            function idxs = getIdxFromFreq(startFreq,stopFreq,freqSamples)
            % Get the frequency sample indexes from start and stop frequencies
                startIdxs = find(startFreq <= freqSamples);
                stopIdxs = find(stopFreq >= freqSamples);
                idxs = (startIdxs(1):stopIdxs(end));
            end
        end
    end
end

function txBandsNormValues = normalizeTxBandObjectiveValues(txBandsObjValues)
    % Normalize tx band objective values
    normColumns = {'SINR','Bandwidth'};
    txBandsNormValues = normalizeColumns(txBandsObjValues,normColumns);
        
    function txBandsNorm = normalizeColumns(txBands,normColumnNames)
        txBandsNorm = txBands;    
        [rows,~] = size(txBands);
        
        for iColumn = 1:length(normColumnNames)
            name = normColumnNames{iColumn};
            normName = ['norm',name];
            maxValue = max([txBands.(name)]);
            for iRow = 1:rows
                txBandsNorm(iRow).(normName) = txBands(iRow).(name)/maxValue;
            end
        end
    end
end

function txBands = getMaxObjectiveBands(normObjective,alpha)
    % Get the idxs that maximize the objective function
    maxObjective = -Inf;
    maxObjectiveRow = 0;
    [rows,~] = size(normObjective);

    for iRow = 1:rows
        objective = alpha*normObjective(iRow).normSINR + (1-alpha)*normObjective(iRow).normBandwidth;
        if objective > maxObjective
            maxObjectiveRow = iRow;
            maxObjective = objective;
        end
    end

    txBands = normObjective(maxObjectiveRow).BandCombo;
end

function notchFrequency = getNotchFrequencyFromTxRegion(mesoBands,txMesoBands,spectrumFrequencies)
    % Get the frequency ranges of the notches based on the selected tx
    % meso-bands
    notchMesoBands = mesoBands;
    notchMesoBands(txMesoBands) = [];

    [~,numNotches] = size(notchMesoBands);
    notchFrequency = zeros(numNotches,2);
    for iNotch = 1:numNotches
        bandStart = notchMesoBands(iNotch).iStart;
        bandStop = notchMesoBands(iNotch).iStop;
        notchFrequency(iNotch,:) = [spectrumFrequencies(bandStart) spectrumFrequencies(bandStop)];
    end
end

end

The following group of helper functions implement the waveform generation algorithms.

function waveform = helperGenerateUnnotchedWaveform(bandwidth,pulseWidth,sampleRate,numPulses)
% Generate an unnotched waveform using the approach described in [1].

% Call helperGenerateNotchedWaveform with empty notch frequencies
notchFrequencies = [];
waveform = helperGenerateNotchedWaveform(bandwidth,pulseWidth,sampleRate,numPulses,notchFrequencies);

end

function waveform = helperGenerateNotchedWaveform(bandwidth,pulseWidth,sampleRate,numPulses,notchFrequencies)
% Generate a notched waveform using the approach described in [1].

% Generate the seed waveform
pcfmWaveform = helperSeedWaveform(bandwidth,pulseWidth,sampleRate,numPulses);

% Calculate the number of samples in the pulse, and get the frequencies in
% the spectrum
numSamples = pulseWidth*sampleRate;

% Generate the unnotched Gaussian spectrum.
unnotchedSpectrum = gausswin(numSamples);

% Notch the spectrum
notchedSpectrum = helperNotchSpectrum(unnotchedSpectrum,sampleRate,notchFrequencies);

% Generate the notched waveform
iterations = 500;
waveform = helperProFm(pcfmWaveform,notchedSpectrum,iterations);

end

function pcfmWaveform = helperSeedWaveform(bandwidth,pulseWidth,sampleRate,numPulses)
% Generate a PCFM waveform with randomized frequencies following the
% process described in [3].

% Get the number of samples in a symbol and the number of symbols in a
% pulse
samplesPerPulse = pulseWidth*sampleRate;
samplesPerSymbol = sampleRate/bandwidth;
symbolsPerPulse = samplesPerPulse/samplesPerSymbol;
samplesPerPulse = samplesPerSymbol*symbolsPerPulse;

% Perform randomized PCFM steps
randomPhases = 2*pi*(rand(symbolsPerPulse+1,numPulses)-0.5);

% Calculate the phase difference
phaseDifference = pcfmCalculatePhaseDifference(randomPhases);

% Generate impulse train
impulseTrain = pcfmGenerateImpulseTrain(phaseDifference,numPulses,samplesPerSymbol,samplesPerPulse);

% Filter impulses
filteredImpulseTrain = pcfmFilterImpulseTrain(impulseTrain,samplesPerSymbol);

% Integrate to get phase
pcfmPhases = randomPhases(1,:) + cumtrapz(filteredImpulseTrain,1);

% Generate samples based on phase
pcfmWaveform = exp(1i*pcfmPhases);

% Get the phase difference of the random phases
function phaseDifference = pcfmCalculatePhaseDifference(randomPhases)
    phaseDifference = diff(randomPhases,[],1);
    greaterThanPi = abs(phaseDifference) > pi;
    phaseDifference(greaterThanPi) = phaseDifference(greaterThanPi) - 2*pi*sign(phaseDifference(greaterThanPi));
end

% Get the impulse train representing frequency
function impulseTrain = pcfmGenerateImpulseTrain(phaseDifference,numPulses,samplesPerSymbol,totalSamples)
    impulseTrain = zeros(totalSamples,numPulses);
    impulseTrain(1:samplesPerSymbol:totalSamples,:) = phaseDifference;
end

% Get the filtered impulse train representing frequency
function filteredImpulseTrain = pcfmFilterImpulseTrain(impulseTrain,samplesPerSymbol)
    raisedCosFilter = rcosdesign(0.25,6,samplesPerSymbol)';

    [rows,cols] = size(impulseTrain);
    filteredImpulseTrain = zeros(rows,cols);
    for iCol = 1:cols
        filteredImpulseTrain(:,iCol) = conv(impulseTrain(:,iCol),raisedCosFilter,'same');
    end
end

end

function proFmWaveform = helperProFm(seedWaveform,desiredSpectrumPower,iterations)
% Perform the PRO-FM technique following the process described in [2].

i = 1;
proFmWaveform = seedWaveform;
desiredSpectrumMagnitude = desiredSpectrumPower.^(1/2);

while i < iterations

    % Get the fft of the current signal
    proFmWaveformSpectrum = fftshift(fft(proFmWaveform));

    % Get the phase of the current signal spectrum
    proFmWaveformSpectrumPhase = angle(proFmWaveformSpectrum);

    % Phase shift the desired spectrum magnitude so that its phase matches
    % the current signal phase
    filteredFrequencySignal = desiredSpectrumMagnitude.*exp(1i*proFmWaveformSpectrumPhase);

    % Get the new signal in the time domain
    filteredTimeSignal = ifft(fftshift(filteredFrequencySignal));

    % Normalize the new signal to be constant amplitude
    proFmWaveform = exp(1i*angle(filteredTimeSignal));

    i = i + 1;
end

end

function notchedSpectrum = helperNotchSpectrum(initialSpectrum,sampleRate,notchFrequencies)
% Insert notches into the given spectrum. Use Tukey Tapers following the
% process described in [3].

% Get the spectrum in dB.
logSpectrum = pow2db(initialSpectrum);

% Get the spectrum frequency samples
spectrumFrequencies = helperSpectrumFrequencies(sampleRate,numel(initialSpectrum));

% Set the notch floor to -85 dB
minValue = -85;

% Get the number of samples in the spectrum
numSamples = length(initialSpectrum);

% taper is B/32, number if samples is 2B, so taper is samples/64
taperLength = ceil(numSamples/64);

% Get the indexes of the notch frequencies in the spectrum
notchIdxs = getNotchStartEndIdx(notchFrequencies,spectrumFrequencies);
[notchNum,~] = size(notchIdxs);

% For each notch, update the spectrum
for iNotch = 1:notchNum
    idxs = notchIdxs(iNotch,:);
    logSpectrum(idxs(1):idxs(end)) = minValue;
    
    % set tukey taper on leading and trailing edges of notch
    fallingTaperStart = idxs(1)-taperLength;
    fallingTaperStartValue = logSpectrum(fallingTaperStart);
    fallingTaper = createTukeyTaper(taperLength,fallingTaperStartValue,minValue);
    logSpectrum(fallingTaperStart:idxs(1)) = fallingTaper;

    risingTaperEnd = idxs(end)+taperLength;
    risingTaperEndValue = logSpectrum(risingTaperEnd);
    risingTaper = createTukeyTaper(taperLength,minValue,risingTaperEndValue);
    logSpectrum(idxs(end):risingTaperEnd) = risingTaper;
end

% Return spectrum power as amplitude
notchedSpectrum = db2pow(logSpectrum);

function taper = createTukeyTaper(numSamples,startValue,endValue)
% Insert a taper where spectral notches occur.
    height = abs(startValue-endValue)/4;
    lowestValue = max([startValue,endValue])-height;
    if startValue > endValue
        flipTaper = false;
    else
        flipTaper = true;
    end

    % create taper
    taper = cos(0.5*pi/numSamples*(0:numSamples));
    
    % stretch taper to height
    taper = taper*height;

    % move taper down to end value
    taper = taper+lowestValue;

    if flipTaper
        taper = flip(taper);
    end
end

% Get the start and end index of the notch frequencies
function notchIdxs = getNotchStartEndIdx(notchFreq,spectrumFrequencies)
    [numNotches,~] = size(notchFreq);
    notchIdxs = zeros(numNotches,2);

    for iNotchFreq = 1:numNotches
        notchStart = notchFreq(iNotchFreq,1);
        notchEnd = notchFreq(iNotchFreq,2);

        aboveNotchStart = find(spectrumFrequencies >= notchStart);
        belowNotchEnd = find(spectrumFrequencies <= notchEnd);
        notchIdxs(iNotchFreq,:) = [aboveNotchStart(1) belowNotchEnd(end)];
    end
end

end

Copyright 2023 The MathWorks, Inc.