[Cyclostationary] explain minimum number of samples for FSM

content/cyclostationary.rst

@@ -285,6 +285,8 @@ Below is a minimal Python implementation of the FSM, which is a frequency-based
     plt.ylabel('Cyclic Frequency [Normalized Hz]')
     plt.show()
 
+Note that due to the way the shift is calculated and rounded to an integer number of samples, it helps to process at least :code:`2 / alpha_resolution` samples at a time.
+
 Let's calculate the SCF for the rectangular BPSK signal we used before, with 20 samples per symbol over a range of cyclic frequencies from 0 to 0.3 using a 0.001 step size:
 
 .. image:: ../_images/scf_freq_smoothing.svg
@@ -903,12 +905,12 @@ OFDM
 
 Cyclostationarity is especially strong in OFDM signals due to OFDM's use of a cyclic prefix (CP), which is where the last several samples of each OFDM symbol is copied and added to the beginning of the OFDM symbol.  This leads to a strong cyclic frequency corresponding to the OFDM symbol length (which is equal to the inverse of the subcarrier spacing, plus CP duration). 
 
-Let's play around with an OFDM signal.  Below is the simulation of an OFDM signal with a CP using 64 subcarriers, 25% CP, and QPSK modulation on each subcarrier.  We'll interpolate by 2x to simulate receiving at a reasonable sample rate, so that means the OFDM symbol length in number of samples will be (64 + (64*0.25)) * 2 = 160 samples.  That means we should get spikes at alphas that are an integer multiple of 1/160, or 0.00625, 0.0125, 0.01875, etc. We will simulate 100k samples which corresponds to 625 OFDM symbols (recall that each OFDM symbol is fairly long).  
+Let's play around with an OFDM signal.  Below is the simulation of an OFDM signal with a CP using 64 subcarriers, 25% CP, and QPSK modulation on each subcarrier.  We'll interpolate by 2x to simulate receiving at a reasonable sample rate, so that means the OFDM symbol length in number of samples will be (64 + (64*0.25)) * 2 = 160 samples.  That means we should get spikes at alphas that are an integer multiple of 1/160, or 0.00625, 0.0125, 0.01875, etc. We will simulate 200k samples which corresponds to 1250 OFDM symbols (recall that each OFDM symbol is fairly long).  
 
 .. code-block:: python
 
     from scipy.signal import resample
-    N = 100000 # number of samples to simulate
+    N = 200000 # number of samples to simulate
     num_subcarriers = 64
     cp_len = num_subcarriers // 4 # length of the cyclic prefix in symbols, in this case 25% of the starting OFDM symbol
     print("CP length in samples", cp_len*2) # remember there is 2x interpolation at the end
@@ -941,19 +943,14 @@ Let's play around with an OFDM signal.  Below is the simulation of an OFDM signa
     n = np.sqrt(np.var(samples) * 10**(-SNR_dB/10) / 2) * (np.random.randn(N) + 1j*np.random.randn(N))
     samples = samples + n
 
-Using the FSM to calculate the SCF at a relatively high cyclic resolution of 0.0001:
-
-.. image:: ../_images/scf_freq_smoothing_ofdm.svg
-   :align: center 
-   :target: ../_images/scf_freq_smoothing_ofdm.svg
-   :alt: SCF of OFDM using the Frequency Smoothing Method (FSM)
-
-Note the horizontal line towards the top, indicating there is a low cyclic frequency.  Zooming into the lower cyclic frequencies, we can clearly see the cyclic frequency corresponding to the OFDM symbol length (alpha = 0.0125).  Not sure why we only get a spike at 2x, and not 1x or 3x or 4x...  Even dropping the resolution by another 10x doesn't show anything else besides the 2x, if anyone knows feel free to use the "Suggest an Edit" link at the bottom of this page.
+Because we expect spikes at 0.00625, 0.0125, 0.01875, we will use a cyclic frequency resolution of 1e-5 so we get an even multiple.  For situations where it's impractical to use such a fine resolution, or the cyclic frequencies are unknown, oversampling can be used (e.g. increasing samples per symbol, in this OFDM example the oversampling factor is 2).  We must also process at least :code:`2 / alpha_resolution` as part of the FSM approach, so 200k samples.  Below are the results, specifically using :code:`alphas = np.arange(0, 0.02, 1e-5)` and max pooling turned on:
 
 .. image:: ../_images/scf_freq_smoothing_ofdm_zoomed_in.svg
    :align: center 
    :target: ../_images/scf_freq_smoothing_ofdm_zoomed_in.svg
-   :alt: SCF of OFDM using the Frequency Smoothing Method (FSM) zoomed into the lower cyclic freqs
+   :alt: SCF of OFDM using the Frequency Smoothing Method (FSM)
+
+Note the three spikes, which would be even more pronounced if we squash RF frequency and plot cyclic frequency in 1D.
 
 External resources on OFDM within the context of CSP:
 

figure-generating-scripts/cyclostationary.py

@@ -10,7 +10,7 @@
 # Simulate Rect BPSK #
 ######################
 
-if True:
+if False:
     N = 100000 # number of samples to simulate
     f_offset = 0.2 # Hz normalized
     sps = 20 # cyclic freq (alpha) will be 1/sps or 0.05 Hz normalized
@@ -153,9 +153,9 @@ def fractional_delay(x, delay):
 ###########################
 
 # Adapted from https://dspillustrations.com/pages/posts/misc/python-ofdm-example.html
-if False:
+if True:
     from scipy.signal import resample
-    N = 100000 # number of samples to simulate
+    N = 200000 # number of samples to simulate
     num_subcarriers = 64
     cp_len = num_subcarriers // 4 # length of the cyclic prefix in symbols, in this case 25% of the starting OFDM symbol
     print("CP length in samples", cp_len*2) # remember there is 2x interpolation at the end
@@ -203,7 +203,7 @@ def fractional_delay(x, delay):
 # Direct CAF #
 ##############
 
-if True:
+if False:
     # CAF only at the correct alpha
     alpha_of_interest = 1/sps # equates to 0.05 Hz
     #alpha_of_interest = 0.08 # INCORRECT ALPHA FOR SAKE OF PLOT
@@ -254,13 +254,13 @@ def fractional_delay(x, delay):
 
 
 # Freq smoothing
-if False:
+if True:
     start_time = time.time()
 
     alphas = np.arange(0, 0.3, 0.001)
-    if False: # For OFDM example
-        #alphas = np.arange(0, 0.5+0.0001, 0.0001) # enable max pooling for this one
-        alphas = np.arange(0, 0.02+0.0001, 0.0001)
+    if True: # For OFDM example
+        #alphas = np.arange(0, 0.5+0.0001, 0.0001) # zoomed out, enable max pooling for this one
+        alphas = np.arange(0, 0.02+0.0001, 0.00001) # zoomed in and increased resolution
     Nw = 256 # window length
     N = len(samples) # signal length
     window = np.hanning(Nw)
@@ -277,15 +277,16 @@ def fractional_delay(x, delay):
     SCF = np.abs(SCF)
 
     # null out alpha= 0, 1, -1 which is just the PSD of the signal, it throws off the dynamic range
-    SCF[np.argmin(np.abs(alphas)), :] = 0 
+    SCF[np.argmin(np.abs(alphas)), :] = 0
+    SCF[np.argmin(np.abs(alphas))+1, :] = 0 # PSD bleeds into 2nd row
     SCF[np.argmin(np.abs(alphas - 1)), :] = 0
     SCF[np.argmin(np.abs(alphas - (-1))), :] = 0
 
     print("Time taken:", time.time() - start_time)
 
     print("SCF shape", SCF.shape)
-    # Max pooling in cyclic domain
-    if False:
+    # Max pooling in cyclic domain, used for OFDM example and possibly others
+    if True:
         import skimage.measure
         SCF = skimage.measure.block_reduce(SCF, block_size=(16, 1), func=np.max) # type: ignore
         print("Shape of SCF:", SCF.shape)
@@ -296,7 +297,7 @@ def fractional_delay(x, delay):
     plt.ylabel('Cyclic Frequency [Normalized Hz]')
     #plt.savefig('../_images/scf_freq_smoothing.svg', bbox_inches='tight')
     #plt.savefig('../_images/scf_freq_smoothing_ofdm.svg', bbox_inches='tight') # for OFDM example
-    #plt.savefig('../_images/scf_freq_smoothing_ofdm_zoomed_in.svg', bbox_inches='tight') # for OFDM example 2
+    plt.savefig('../_images/scf_freq_smoothing_ofdm_zoomed_in.svg', bbox_inches='tight') # for OFDM example 2
     #plt.savefig('../_images/scf_freq_smoothing_pulse_shaped_bpsk.svg', bbox_inches='tight')
     #plt.savefig('../_images/scf_freq_smoothing_pulse_shaped_bpsk2.svg', bbox_inches='tight')
     #plt.savefig('../_images/scf_freq_smoothing_pulse_shaped_bpsk3.svg', bbox_inches='tight')