{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#A worksheet to analyze frequency content of a signal from sampled data\n",
    "#using the discrete cosine transform (DCT).\n",
    "#Parameters tlow and thigh control plot range for time domain plot.\n",
    "#Parameters flow and fhigh control frequency range for DCT plot.\n",
    "\n",
    "#First load in the appropriate data set, here in the form of (time,signal) pairs\n",
    "#from the supplied files.\n",
    "\n",
    "#load('gong.sage'); samprate = 16000.0; #sampled a 16 kHz\n",
    "#load('sunspot_data.sage'); samprate = 1.0; #sampled at 1 per year\n",
    "load('spring_mass_data_fourier.sage'); samprate = 50.0; #sampled at 50 Hz."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The data is in the array \"udat0\".\n",
    "#Number of data points \"N\" and time interval \"T\":\n",
    "N = len(udat0);\n",
    "T = N/samprate;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Select a time interval for plot and plot.\n",
    "tlow = 0.0; thigh = T;\n",
    "jlow = floor(N*tlow/T); jhigh = min(ceil(N*thigh/T),N);\n",
    "J2 = list(range(jlow,jhigh));\n",
    "pdat = [udat0[j] for j in J2];\n",
    "scatter_plot(pdat,facecolor='red',markersize=10,axes_labels=['Time','Signal Magnitude']).show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Compute the DCT of the signal, bootstrapping off of the fft command.\n",
    "\n",
    "J = list(range(N)); #list 0 to N-1\n",
    "A = [RR(udat0[i][1]) for i in J]; #Pick off signal components\n",
    "A2 = [A[i] for i in J]; #Make a copy\n",
    "A.reverse(); #Reverse A\n",
    "A2.extend(A); #and append to A2\n",
    "J2 = list(range(2*N));\n",
    "udat = IndexedSequence(A2,J2); #Data structure for FFT command\n",
    "C = udat.fft(); #Compute FFT\n",
    "fC = C.list(); #Pick off components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Modify fC to obtain DCT of data.\n",
    "p2 = n(pi/2);\n",
    "W = [exp(-p2*I*j) for j in J];\n",
    "DC = [real(W[j]*fC[j])/sqrt(2.0*n(N)) for j in J];\n",
    "DC[0] = DC[0]/sqrt(2.0); #Rescale DC[0] special."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Prepare to plot DCT coefficient magnitudes.\n",
    "J = list(range(0,N));\n",
    "DCmag = [abs(DC[j]) for j in J];\n",
    "freqs = [0.5*j*samprate/n(N) for j in J];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Select frequency range for plot, minimum frequency 0, maximum samprate/2\n",
    "flow = 1.0; fhigh = 25.0;\n",
    "jlow = floor(2*flow*N/samprate); jhigh = min(ceil(2*fhigh*N/samprate),N);\n",
    "J2 = list(range(jlow,jhigh));\n",
    "pdat = [[freqs[j],DCmag[j]] for j in J2];\n",
    "scatter_plot(pdat,facecolor='red',markersize=10,axes_labels=['Frequency (Hz)','DCT Magnitude']).show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.2",
   "language": "sage",
   "name": "sagemath"
  },
  "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}