{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c3f7f94a-e06d-425e-b68f-1116ce2263d2",
   "metadata": {},
   "source": [
    "# Complex Ambiguity Function (CAF)\n",
    "The complex ambiguity function is defined as:\n",
    "\n",
    "${\\displaystyle \\chi (\\tau ,f)=\\int _{-\\infty }^{\\infty }s(t)s^{*}(t-\\tau )e^{i2\\pi ft}\\,dt}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4feeaf64-f0a0-4407-8840-75f772c69900",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import cm\n",
    "\n",
    "import numpy as np\n",
    "from scipy import signal\n",
    "from sk_dsp_comm import sigsys as ss\n",
    "from sk_dsp_comm import fir_design_helper as fir_h\n",
    "\n",
    "%config InlineBackend.figure_formats=['svg'] # SVG inline viewing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5702bef0-b3fd-49ba-8323-6c375a2be665",
   "metadata": {},
   "source": [
    "### FFT CAF using Transform Domain "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ed6dd8f-68fc-433d-b68e-3cdbc39d9239",
   "metadata": {},
   "source": [
    "#### Use a Lowpass Filtered White Gaussian Noise Signal\n",
    "Here we design a lowpass filter to create a waveform with correlation. We snip out samples 900 to 1200 to form a correlation reference waveform. Stream in 2048 samples and use an FFT half-length of 1024. Display the CAF magnitude as an image and as a 3D surface. We also include a 10 Hz shift of the input to verify that the CAF correlation peak is centered at 10 Hz.\n",
    "\n",
    "**Note** you have to make sure `Nfft2 > len(b_lpf)` so the transform domain correlation will work properly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4637b42d-ff71-4fbf-8cf2-5708632b4f15",
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = 1000\n",
    "fd = 10\n",
    "b_lpf = fir_h.fir_remez_lpf(100,150,0.1,80,fs)\n",
    "print('N_b_lpf = %d' % len(b_lpf))\n",
    "xp = signal.lfilter(b_lpf,1,np.random.randn(2048))\n",
    "xp_ref = xp[900:1200] # create a reference waveform from xp\n",
    "n = np.arange(len(xp))\n",
    "yp = xp * np.exp(1j*2*np.pi*fd/fs*n)\n",
    "# Use a non-power of 2 FFT to exact slice spacings (FFTW is somewhat slower)\n",
    "y_caf_stream,faxis,taxis = ss.fft_caf(yp,xp_ref,n_fft2=1000,n_slice2=80,bs=0.5,fs = 1000.0)\n",
    "plt.figure() # figsize=(6,6))\n",
    "plt.imshow(abs(y_caf_stream[:,1150:1250]), extent=[taxis[1150], taxis[1250], faxis[-1], faxis[0]], aspect='auto')\n",
    "plt.ylabel(r'Frequency (Hz)')\n",
    "plt.xlabel(r'Time (seconds)')\n",
    "plt.title(r'CAF Magnitude')\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45a76ffc-b12b-447c-ad58-eb517c34d1ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a 3D figure with axes with a specified size\n",
    "fig = plt.figure(figsize=(8, 8))  # Adjust the width and height as needed\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "# fig, ax = subplots(subplot_kw={\"projection\": \"3d\"})\n",
    "tgrid, fgrid = np.meshgrid(taxis[1150:1250], faxis)\n",
    "surf = ax.plot_surface(tgrid, fgrid, abs(y_caf_stream[:,1150:1250]), cmap=cm.jet,\n",
    "                       linewidth=0, antialiased=False)\n",
    "plt.ylabel(r'Frequency (Hz)')\n",
    "plt.xlabel(r'Time (seconds)')\n",
    "plt.title(r'CAF Magnitude')\n",
    "# Change the view angle\n",
    "ax.view_init(elev=30, azim=135)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09f45ed5-d6d4-4e33-9747-da247ba8fea7",
   "metadata": {},
   "source": [
    "#### Use a PN Sequence\n",
    "Here we use a PN6 for correlation with a 10,000 bit, lowpass filtered, NRZ waveform. The bit rate is  1 bps. The correlation waveform and input waveform each have 10 samples per bit. The reference waveform is pure NRZ, with the filtering only present on the input waveform. The FFT half-length is 1024 to increase the frequency resolution to $f_s/2048 = 10/2048 = 0.00488$ Hz."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52438da7-b94b-4e55-b3b9-827411dcd242",
   "metadata": {},
   "outputs": [],
   "source": [
    "b_nrz = fir_h.fir_remez_lpf(1.25,1.8,0.1,50,10)\n",
    "pn_data = np.int16(ss.pn_gen(10000,6))\n",
    "x_nrz, b2 = ss.nrz_bits2(pn_data,10,'rect')\n",
    "# Reduce the Sidelobe level of the NRZ waveform\n",
    "y_nrz = signal.lfilter(b_nrz,1,x_nrz)\n",
    "Px_nrz, fx_nrz = ss.psd(y_nrz,2**10,10)\n",
    "plt.plot(fx_nrz,10*np.log10(Px_nrz))\n",
    "plt.title(r'Filtered NRZ Input Signal')\n",
    "plt.ylabel(r'PSD (dB)')\n",
    "plt.xlabel(r'Frequency (Hz)')\n",
    "plt.grid();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8bbff39d-0ab4-41b6-920e-8422f91281f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "pn_period = np.int16(ss.m_seq(6))\n",
    "h_ref, b1 = ss.nrz_bits2(pn_period,10,'rect')\n",
    "h_ref2 = signal.lfilter(np.ones(5)/5,1,h_ref)\n",
    "pn_data = np.int16(ss.pn_gen(300,6))\n",
    "x_nrz, b2 = ss.nrz_bits2(pn_data,10,'rect')\n",
    "y_nrz = signal.lfilter(b_nrz,1,x_nrz)\n",
    "n = np.arange(len(y_nrz))\n",
    "y_nrz_s = y_nrz * np.exp(1j*2*np.pi*(-2)/fs*n)\n",
    "y_caf_stream2,faxis2,taxis2 = ss.fft_caf(y_nrz_s,h_ref2,n_fft2=1000,n_slice2=15,bs=1,fs = 1000.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8018a6e-c374-4d57-9950-03fea9a7a581",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure() # figsize=(6,6))\n",
    "plt.imshow(abs(y_caf_stream2[:,500:800]), extent=[taxis2[500], taxis2[800], faxis2[-1], faxis2[0]], aspect='auto')\n",
    "plt.ylabel(r'Frequency (Hz)')\n",
    "plt.xlabel(r'Time (seconds)')\n",
    "plt.title(r'CAF Magnitude')\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f222474-7682-48d2-819a-f530a161c9d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a figure and axes with a specified size\n",
    "fig = plt.figure(figsize=(8, 8))  # Adjust the width and height as needed\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "# fig, ax = subplots(subplot_kw={\"projection\": \"3d\"})\n",
    "tgrid, fgrid = np.meshgrid(taxis2[500:800], faxis2)\n",
    "surf = ax.plot_surface(tgrid, fgrid, abs(y_caf_stream2[:,500:800]), cmap=cm.jet,\n",
    "                       linewidth=0, antialiased=False)\n",
    "plt.ylabel(r'Frequency (Hz)')\n",
    "plt.xlabel(r'Time (seconds)')\n",
    "plt.title(r'CAF Magnitude')\n",
    "# Change the view angle\n",
    "ax.view_init(elev=40, azim=145)\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}