{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Research example\n",
    "This notebook preserves a 2019 exploratory method by Andrew Dyck.\n",
    "Read README.md before running it. Full scientific validation is still needed.\n",
    "Outputs below have been cleared and do not represent a new execution.\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "from pathlib import Path\n",
    "Path(\"generated\").mkdir(exist_ok=True)\n"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# FRB Pipeline Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "import numpy as np\n",
    "from matplotlib.pyplot import *\n",
    "import gc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def axis_labels(x, y, z):\n",
    "    xlab = xlabel(x)\n",
    "    ylab = ylabel(y)\n",
    "    titles = title(z)\n",
    "    legends = legend\n",
    "    return xlab, ylab, titles, legends\n",
    "\n",
    "def centfinder(x, threshold, length):\n",
    "    intensity = x\n",
    "    peaksx = [] #This is the values that the peaks are centred around\n",
    "    peaksy = []\n",
    "    peak_dex = (list(intensity))\n",
    "    inversep_dex=intensity[::-1]\n",
    "    for i in range(len(intensity)-1):\n",
    "        if intensity[i] > threshold and intensity[i-2]<intensity[i] and intensity[i-1] <intensity[i] and intensity[i+1] < intensity[i] and intensity[i+2] < intensity[i]:\n",
    "            peaksx.append(length[i])\n",
    "            peaksy.append(intensity[i])\n",
    "                \n",
    "    #Calculate the FWHM to compute the size of the peak that we will care about, and then the centroid\n",
    "    cen = []\n",
    "    for i in range(len(peaksy)):\n",
    "        halfmax = peaksy[i]/2\n",
    "        peak_index = peak_dex.index(peaksy[i])\n",
    "        pind=(list(inversep_dex)).index(peaksy[i])\n",
    "        \n",
    "        for j in range(len((intensity[0:pind]))):\n",
    "            xmin = (inversep_dex[pind:-1][j] - halfmax)\n",
    "            xmins = peak_index-j+1\n",
    "            if xmin < 0 or inversep_dex[pind:-1][j] > peaksy[i]:\n",
    "                break\n",
    "\n",
    "        for j in range(len((intensity[0:peak_index] - halfmax))):\n",
    "            xmax = (intensity[peak_index:-1][j] - halfmax)\n",
    "            xmaxs = peak_index+j-1\n",
    "            if xmax < 0 or intensity[peak_index:-1][j] > peaksy[i]:\n",
    "                break\n",
    "        \n",
    "        x_range = length[xmins:xmaxs]\n",
    "        y_range = x[xmins:xmaxs]\n",
    "    return x_range, yrange\n",
    "\n",
    "\n",
    "def peakfinder(x, threshold, array_length):\n",
    "    \"\"\"This function finds all maxima above a certain threshold, with the condition that the two points before\n",
    "    and after the maxima must have values below it.\"\"\"\n",
    "    pixels = array_length\n",
    "    intensity = x\n",
    "    peaksx = [] #This is the values that the peaks are centred around\n",
    "    peaksy = []\n",
    "    peak_dex = (list(intensity))\n",
    "    inversep_dex=intensity[::-1]\n",
    "    for i in range(len(intensity)-1):\n",
    "        if intensity[i] > threshold and intensity[i-2]<intensity[i] and intensity[i-1] <intensity[i] and intensity[i+1] < intensity[i] and intensity[i+2] < intensity[i]:\n",
    "            peaksx.append(pixels[i])\n",
    "            peaksy.append(intensity[i])\n",
    "    return peaksx, peaksy\n",
    "\n",
    "\n",
    "def centfinder_cheat(x, threshold):\n",
    "    \"\"\"This works well for one spike, not multiple. What is does is it find the peak values and then\n",
    "    calculates the centroid for all the peaks\"\"\"\n",
    "    pixels = velocities\n",
    "    intensity = x\n",
    "    peaksx = [] #This is the values that the peaks are centred around\n",
    "    peaksy = [np.max(intensity)]\n",
    "    peak_dex = (list(intensity))\n",
    "    inversep_dex=intensity[::-1]\n",
    "                \n",
    "    #Calculate the FWHM to compute the size of the peak that we will care about, and then the centroid\n",
    "    cen = []\n",
    "    for i in range(len(peaksy)):\n",
    "        halfmax = peaksy[i]/2\n",
    "        peak_index = peak_dex.index(peaksy[i])\n",
    "        pind=(list(inversep_dex)).index(peaksy[i])\n",
    "        \n",
    "        for i in range(len((intensity[0:pind]))):\n",
    "            xmin = (inversep_dex[pind:-1][i] - halfmax)\n",
    "            xmins = peak_index-i+1\n",
    "            if xmin < 0:\n",
    "                break\n",
    "\n",
    "        for i in range(len((intensity[0:peak_index] - halfmax))):\n",
    "            xmax = (intensity[peak_index:-1][i] - halfmax)\n",
    "            xmaxs = peak_index+i-1\n",
    "            if xmax < 0:\n",
    "                break\n",
    "        \n",
    "        x_range = velocities[xmins:xmaxs]\n",
    "        y_range = x[xmins:xmaxs]\n",
    "\n",
    "        frequency = np.sum(x_range*y_range)/ np.sum(y_range)\n",
    "        frequency\n",
    "        cen.append(frequency)\n",
    "    return cen"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following steps will simulate the anologue process that occurs before the data is downsampled. Any pulsar will work for this, however, the pulsar being simulated will be made using the paramters below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Define the parameters for the noise that will be mixed with the pulse\n",
    "Cf = 1400e6 #This is the centre frequency in MHz, and where we will be viewing from\n",
    "Bw = 10e6 #This is the bandwidth that we are observing with our telescope\n",
    "sample_rate = 12e9 # Frequency running simulation at.\n",
    "period = .005 #seconds "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With a period of 5ms, we will use a time frame with this length to simulate one pulse. The bins in this simulation is determined by the sampling rate, which can be changed above. The current value is at such that 12GB ram computerss will just be able to simulate. The higher value the sampling rate is given, the more accurate the simulation becomes, but at the expense of what the computer is able to process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Define the timescale in which we will be observing\n",
    "t = np.linspace(-.0025,.0025, int(period*sample_rate)) #Decreasing the period for this will decrease the samples obtained"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To create the pulsar, first noise is simulated using a normal distribution of values. This distribution is made such that the standard deviation is 1, and is centred around 0. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Create the Gaussian noise profile for the pulsar\n",
    "length = len(t)     #This defines the number of timesamples we have, ie the number of samples collected\n",
    "mu = 0              #This indicates an offset of 0, meaning the signal is perfetly centred about 0\n",
    "sigma = 1           #This idicates an std of 1\n",
    "gnoise = np.random.normal(mu, sigma, length) #This is the GAussian noise for the signal\n",
    "pw = .0004 #seconds\n",
    "print(length)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This noise is then multiplied to a gaussian that has a standard deviation equal to that of the desired pulse width.The gaussian was made using the equation \\begin{equation} e^{\\frac{-1}{2}(\\frac{t}{pw})^2} \\end{equation}."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gnoise = np.exp(.5*-(t/pw)**2)*gnoise #Define the pulse of the signal. The gaussian profile takes an input of the "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observe now that the resulting signal looks like a Gausian distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now plot the signal, and observe the pulse shape. To see the pulse without destroying the computer, \n",
    "#plot every 10000th point\n",
    "figure()\n",
    "plot(np.linspace(0,.005, int(period*sample_rate))[0::10000],gnoise[0::10000], ls='none', marker='.')\n",
    "axis_labels('time(s)', 'value(unitless)', 'Gaussian noise pulse')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, to create the pulsar, we take the Fourier transform of the Gaussian, so that we may eliminate all but one frequency of the Gaussian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now take the fourier transform of the pulse to get it in frequency space\n",
    "fourier_pulse = np.fft.fft(gnoise)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The x-axis for the resulting graph is caluclated here. This tells what the frequency step for the broadband noise of the Gaussian is. This is done by taking the fftfreq of the length of the Gaussian, and having it step by the sample rate in units of MHz."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The fourier transformed data should look like broadband noise, so define the axes properly before plotting\n",
    "broad_freq = np.fft.fftfreq(length, 1/(sample_rate/1e6)) #This gives the signal in MHz, (specifically from the 1e6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot of the Fourier transform of the Gaussian wave packet is a broadband noise plot as expected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "figure()\n",
    "plot(broad_freq[::1000], np.abs(fourier_pulse[::1000])) #Square the abs value to get better view\n",
    "axis_labels(\"frequency(MHz)\", 'amplitude(unitless)', 'Fourier transformed  gaussian pulse')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, all frequencies but one are reduced to zero through filtering. To filter, the highpass frequency, or the highest frequency that we wish to observe with our bandwidth, is used to find the point in our data that is set as an index....."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now we modify the noise, so that there are only +- of the same frequency in it of it. This is the pulsar!\n",
    "#Everything outside the two reference frequencies is turned to zero\n",
    "highpass_freq = (Cf-Bw/2.0) #This sets the centre frequency to 0, with a bandwidth deifned above\n",
    "highpass_index1 = int(highpass_freq/(sample_rate/2)*length/2) #This finds the point in the array with the value\n",
    "highpass_index2 = length - highpass_index1 #Same but for different point\n",
    "fourier_pulse[:highpass_index1] = 0 #This is turning everything outside the two points to zero\n",
    "fourier_pulse[highpass_index2:] = 0 #Same as above\n",
    "\n",
    "#Everything outside of the two frequencies is turned to 0\n",
    "nyquest_index = int(length/2) #This value is determined by the nyquist sampling theorem(helps figure out where the chinkc are)\n",
    "#The following is the same as above, but now for outside the frequency values\n",
    "lowpass_freq = Cf+Bw/2.0 \n",
    "lowpass_index1 = int(lowpass_freq/(sample_rate/2)*length/2)\n",
    "lowpass_index2 = length - lowpass_index1\n",
    "fourier_pulse[lowpass_index1:nyquest_index] = 0\n",
    "fourier_pulse[nyquest_index:lowpass_index2] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "figure()\n",
    "plot(broad_freq[::1000], np.abs(fourier_pulse[::1000]))\n",
    "axis_labels(\"frequency(MHz)\", 'amplitude(unitless)', 'Newly created pulsar')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that all but the absolute value of one frequency has been given a value of zero, we have created our pulsar. The pulsar is now changed back into the time domain using an inverse fourier transform, and still looks like a Gaussian wave packet as desired. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Now turn the pulsar back into time space from frequency space\n",
    "blimited_pulse = np.fft.ifft(fourier_pulse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#If this pulse is plotted again, it sill looks like the same as before\n",
    "figure()\n",
    "plot(np.linspace(0,.005, int(period*sample_rate))[0::10000], blimited_pulse[::10000])\n",
    "axis_labels('time(s)', 'value(unitless)', 'Pulsar wave packet')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To simulate what the computer recieves after the telescope detects the signal, we need to create a mixing signal which will be used to downsample the data. This mixing signal is created using the complex Euler equation:\n",
    "$e^{-2\\pi j *\\frac{Cf}{samp}*len}$, where Cf is defined as above, samp is the sample rate, and len is the array with size length, with step size of 1. It is also assumed that the pulse currently is centred around the same point, x=0, which eliminates the \\omega that is in the original equation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now need to create the mixing signal:\n",
    "mixing_signal = np.exp(-2.00j*np.pi* Cf/sample_rate *np.arange(length)) #With this we create the mixing signal using the euler equation for the signal\n",
    "\n",
    "#We assume the pulse is always centred around the same point, so the only part of the signal equaiton needed is \n",
    "#omega*t. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The pulse is then multiplied to the mixing signal. The fourier transform of this new signal is then taken, and a lowpass filter is applied to it. This is to create an IQ data set for the bandwidth. This filter was made using the nyquist sampling theorem to prevent oversampling with the bandwidth specified. The inverse fourier transform of this signal is then taken, and it is downsampled by taking only samples that have spcings of the $\\frac{sample\\_ rate}{bandwidth}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "#Mix and sample down the pusar signal so that we can have an I Q data set for the bandwidth\n",
    "mixedd_down_s = blimited_pulse*mixing_signal\n",
    "N_cutoff = int(Bw/2/sample_rate*length) #This is from the lowpass fileter definition\n",
    "fmixed = np.fft.fft(mixedd_down_s) #This creates the fourier transform of the mixed signal\n",
    "fmixed[N_cutoff:-N_cutoff] = 0.0 # this makes it so that the imaginary magnitudes are comparable to the reals. This is\n",
    "                                    #also our filter\n",
    "filtered_mixed_downs_s = np.fft.ifft(fmixed) #This brings the fourier space pulse back into a packet\n",
    "#Downsample so that is...\n",
    "downsampled_filtered_mixed_down_s = filtered_mixed_downs_s[::int(sample_rate/(Bw))]\n",
    "#complex sampled at 10MHz instead of 12GHz\n",
    "#This gives us 50000 indicies instead of 60000000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To open up more RAM in the computer, variables that will no longer be used are forgotten"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%reset_selective -f 't'\n",
    "%reset_selective -f gnoise \n",
    "%reset_selective -f fourier_pulse\n",
    "%reset_selective -f broad_freq\n",
    "%reset_selective -f blimited_pulse\n",
    "%reset_selective -f mixing_signal\n",
    "%reset_selective -f mixedd_down_s\n",
    "%reset_selective -f N_cutoff\n",
    "%reset_selective -f fmixed\n",
    "%reset_selective -f filtered_mixed_downs_s\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.pyplot import *\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The pulse is now downsampled, and realistically, this is where the simulation can begin. The previous lines of code were used to demponstrate how the anologue to digital converter chain works. Now, the new downsampled signal will be plotted below. Again, this pulsar signal still looks like a wave packet. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#The simulation can be run from this point\n",
    "\n",
    "\n",
    "# Use gaussian noise, filter with same time width.\n",
    "#Now plot the signal, and observe that the pulse shape is the same, but now using smaller increments between samples\n",
    "figure()\n",
    "plot(np.linspace(0,.005, len(downsampled_filtered_mixed_down_s))[0:], downsampled_filtered_mixed_down_s[:].real)\n",
    "plot(np.linspace(0,.005, len(downsampled_filtered_mixed_down_s))[0:], downsampled_filtered_mixed_down_s[:].imag)\n",
    "axis_labels('Time(s)', 'Amplitude', 'Downsampled Pulsar')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Right now, there is a single pulse, so this is more of an FRB than a pulsar. In order to make this a repeating pulse, use numpy to repeat it. Below, this was done 30 times to get a nice sample of pulses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#We repeat the pulse to get about 0.15 seconds of data\n",
    "downsampledd_filtered_mixed_down_s_repeated = np.tile(downsampled_filtered_mixed_down_s, 30)\n",
    "time_length = .005*30"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now finally, the pulsar is dispersed, which is what the telescope would show the pulsar as. The number of time samples of the new repeated pulse is found, and this is used to find the frequency steps for the pulse. A Dispersion measure (DM) is arbitrarilly chosen. Then the fourier transform of the delay time, given by the equation\n",
    "\\begin{equation}4.15e6 *DM*(\\frac{1}{\\nu_1^2}-\\frac{1}{\\nu_2^2}) \\end{equation}. This equation assumes a frequency given in MHz, however here, it is given in Hz. Therefore, the equation must be divided by 1e9, which is how the constant in the code below is obtained. The fourier transform of this is then taken and convolved to the fourier transform of the pulse to disperse it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#This is where the pulsar is dispersed\n",
    "#We do a convolution in fourier space, which must be long wrt the dispersion delay time\n",
    "dis_length = len(downsampledd_filtered_mixed_down_s_repeated) #This gives us the number of time samples of the \n",
    "                                                             #downsampled repeated pulse\n",
    "freq_step = np.linspace(-Bw/2, Bw/2, dis_length) #This is the new frequency step for the pulse\n",
    "DM = 600 #This is the dispersion measure and it is completely arbitrary\n",
    "\n",
    "#Create the signal to be convolved with the data. This signal is what disperses the pulsar\n",
    "H = np.exp(2j*np.pi*4.15e15*DM*freq_step**2/((Cf+freq_step)*Cf**2)) #Where the part that looks similar to the delay \n",
    "#time is the 'time' portion of the equation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a precaution, the delay time is checked to ensure that it is less than the dispersion time of the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Check how long the dispersion is\n",
    "delay_time = 4.15e15 * (1/(Cf-Bw/2)**2 - 1/(Cf+Bw/2)**2) * DM\n",
    "\n",
    "#The delay time must be significantly less than the dispersion of the data\n",
    "print(delay_time) #This time must be less thah H"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The pulse signal and the dispersion function are convolved together, and the outcome is the dispersed pulse train"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dispersed_downsampled_filtered_mixed_sig = np.fft.ifft(np.fft.fft(downsampledd_filtered_mixed_down_s_repeated)*np.fft.fftshift(H))\n",
    "#The pulse signal and the dispersion of the pulse were convolved to give a dispersed pulse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The pulse has successfully been dispersed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#Plotting the convolved signal...\n",
    "figure()\n",
    "plot(np.linspace(0,.15, dis_length)[::50], dispersed_downsampled_filtered_mixed_sig[::50].real)\n",
    "plot(np.linspace(0,.15, dis_length)[::50], dispersed_downsampled_filtered_mixed_sig[::50].imag)\n",
    "axis_labels('time(s)', 'Amplitude', 'Dispersed downsampled mixed pulsar (DM=600) ')\n",
    "savefig('generated/despirsed_pulse.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To simulate a real pulse, noise is added to the pulse train. This is white noise generated by choosing a random selection of values within a amplitude range of the pulse,and adding imaginary and real noise to the pulses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Noise is added to the pulse to simulate the RFI environment\n",
    "noise1 = np.random.normal(0,.02,len(dispersed_downsampled_filtered_mixed_sig),)\n",
    "noise2 =  np.random.normal(0,.02,len(dispersed_downsampled_filtered_mixed_sig),)*1j\n",
    "dispersed_downsampled_filtered_mixed_sig = dispersed_downsampled_filtered_mixed_sig+noise1+noise2\n",
    "all_noise = noise1+noise2\n",
    "\n",
    "pre_observing_time = 2000\n",
    "noise_and_signal = np.concatenate([dispersed_downsampled_filtered_mixed_sig,all_noise])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The new pulse with noise is shown below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "figure()\n",
    "plot(np.linspace(0,.15, dis_length)[::50], dispersed_downsampled_filtered_mixed_sig[::50].real)\n",
    "plot(np.linspace(0,.15, dis_length)[::50], dispersed_downsampled_filtered_mixed_sig[::50].imag)\n",
    "axis_labels('time(s)', 'Amplitude', 'Dispersed downsampled mixed pulsar (DM=600) ')\n",
    "savefig('generated/despirsed_pulse_with_noise.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this point, the pulse can be saved to a file and used as a simulated pulse inside GNUradio. The "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dispersed_downsampled_filtered_mixed_sig.dtype\n",
    "rmax = dispersed_downsampled_filtered_mixed_sig.real.max()\n",
    "print(rmax)\n",
    "immax = dispersed_downsampled_filtered_mixed_sig.imag.max()\n",
    "print(immax)\n",
    "int_max = 2**15-2\n",
    "print(int_max)\n",
    "dmax = max(rmax,immax)\n",
    "print(dmax)\n",
    "dispersed_downsampled_filtered_mixed_signal = dispersed_downsampled_filtered_mixed_sig*int_max/dmax\n",
    "dispersed_downsampled_filtered_mixed_down_s_real = dispersed_downsampled_filtered_mixed_signal.real.astype(np.int16)\n",
    "dispersed_downsampled_filtered_mixed_down_s_imag = dispersed_downsampled_filtered_mixed_signal.imag.astype(np.int16)\n",
    "out = np.zeros(2*len(dispersed_downsampled_filtered_mixed_down_s_real), dtype=np.int16)\n",
    "out[::2] = dispersed_downsampled_filtered_mixed_down_s_real\n",
    "out[1::2] = dispersed_downsampled_filtered_mixed_down_s_imag\n",
    "out.tofile('generated/pulse_sim_10mhz_int16_5ms_period_60dm_1400MHz_center_150ms_long.bin')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now this is how this is how the pulse looks when it enters into the computer. The task is to now dedisperse it. We want an integration time that is dependent on the bandwidth. Therefore the inverse of the bandwidth is taken as the time that comes from the bandwidth. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(400/6000)\n",
    "print(3000000/6000)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#We want a 10Msps, since this is the bandwidth\n",
    "times = 1/10e6\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The number of desired frequency channels is set. This is an arbitrary choice. The integration size determins the signal to noise ratio of the data,and so an integration number of .2ms is chosen, since it is much smaller than the pulse width. Therefore, the resulting number of integrations done on the pulsar is 5. Finally, the number of timesteps for the pulse is chosen. This will dictate how many pulses will be able to be seen.\n",
    "\n",
    "With these parameters, the FFT of the data per timestep is done, and arranged into an array as frequency by time.\n",
    "In order to get this into a more conventional array, its transpose is taken, so that it is arranged frequency by time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "#To check that this number of frequencies, integrations, and sampling rate work\n",
    "#integrate 5 samples\n",
    "integration_time = 2e-4 #seconds, this is the 1/.2ms integration\n",
    "step = 400 #This comes from the number of frequency channels desired. This is an educated, but arbitrary choice\n",
    "integration_size = int(integration_time/(times)/step) #This comes from the 0.2ms integration that we desire. It is .2 bc .1 was too small\n",
    "n_time = int(len(dispersed_downsampled_filtered_mixed_sig)/step/integration_size) #This finds the number of\n",
    "                                                            #timesteps\n",
    "spects = []\n",
    "rec_spect = np.zeros(((n_time), step), dtype=complex) #This the is the recieved 3d plot of frequency vs time\n",
    "for i in range(n_time): #This is the number of 'rows'\n",
    "    for j in range(integration_size):\n",
    "        spect = np.fft.fft(dispersed_downsampled_filtered_mixed_sig[(integration_size*i+j)*step : (integration_size*i+j+1)*step], 400)\n",
    "        spects.append(spect)\n",
    "        #This is an easy way to get the row size that we desire\n",
    "        rec_spect[i] += spect*spect.conjugate() #The first timestep of the frequencies needs to be multiplied by its conjugate\n",
    "        #The spetra are added together because we are 'folding' the amplitudes together. If we don't, then there appears \n",
    "        #a jump in the data.\n",
    "\n",
    "        \n",
    "        \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "200/8000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "400*5*1.42*5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#figure()\n",
    "#plot(out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For bench testing purposes, the pulse that has now gone through an FFT and been integrated will be compared to the integration and FFT blocks in GNUradio. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rmax_FFT = rec_spect.real.max()\n",
    "print(rmax)\n",
    "immax_FFT = rec_spect.imag.max()\n",
    "print(immax)\n",
    "int_max = 2**15-2\n",
    "print(int_max)\n",
    "dmax_FFT = max(rmax_FFT,immax_FFT)\n",
    "recieved_spect = (np.fft.fftshift(rec_spect*int_max/dmax_FFT))\n",
    "rec_spect_real = recieved_spect.real.astype(np.int16)\n",
    "rec_spect_imag = recieved_spect.imag.astype(np.int16)\n",
    "out_spec = np.zeros((len(rec_spect_real),step), dtype=np.int16)\n",
    "for i in range(step):\n",
    "    out_spec[:750,i] = rec_spect_real[:,i]\n",
    "out_spec.tofile('generated/pulse_sim_10mhz_int16_5ms_period_60dm_1400MHz_center_150ms_long_integrated_FFT.bin')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then the x axis of the plot is created, where first the number of frequency channels is found, then the centre frequency of that array is set at 1400MHz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_frequencies = np.fft.fftfreq(step, 1/(10e6/1e6)) #This is the number of frequency channles.The first parameter \n",
    "                                #should be the number of chanels desired, and the second is the sample rate in MHz.\n",
    "                                #This will lable the the x axis\n",
    "pulse_freqs = np.fft.fftshift(n_frequencies)+Cf/1e6 #This will place the centre frequency of the plot at 1400MHz\n",
    "dispersed = np.flip(np.transpose(np.fft.fftshift((rec_spect))), 0)#The array of the data must be flipped to match \n",
    "                                    #the units in the plot,and it must be transposed so that time is on the x axis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The 3d graph of the dispersed pulsar is then plotted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#This is the graph of the pulsar, where we can see that it is dispersed\n",
    "figure()\n",
    "imshow(abs(dispersed), extent=[np.amin(np.linspace(0,.15, dis_length)),\\\n",
    "                            np.amax(np.linspace(0,.15, dis_length)),\\\n",
    "                            np.amin(pulse_freqs),\\\n",
    "                            np.amax(pulse_freqs)],\\\n",
    "                              aspect='auto')\n",
    "axis_labels('time(s)','Frequency (MHz)',  'Dispersed Pulsar (DM=600)')\n",
    "colorbar(label = 'power')\n",
    "#savefig(\"generated/Dispersed_pulsar_DM=600.png\")\n",
    "#Congrats!, we now have a dispersed pulsar\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dispersed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the pulsar is attempted to be dedispersed. Assume that the DM, nor the period of the pulsar is known. An array of possible DM's is created, along with a timestep array of when the pulses arrive. Then using the time delay equation used to disperse the pulsar, the time changes for every frequency are calculated. This is repeated for each trial dispersion measure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "figure()\n",
    "plot((out_spec)[3])\n",
    "\n",
    "axis_labels('time(s)','amplitude','Dipersed pulse at one frequency over time')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rec_spect"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DMs = np.linspace(0,400,401) #Assume that we do not know the DM for the pulsar. We must then Guess a range of DM\n",
    "\n",
    "time_change = np.zeros((len(DMs), step)) #This is the time after the first frequency which the second frequncy arrives\n",
    "for i in range(len(DMs)):\n",
    "    change = 4.148e3*DMs[i]*( 1/(np.min(pulse_freqs))**2 - 1/( pulse_freqs )**2) #This is the definition for the time that the second freq arrives\n",
    "    time_change[i,:] = change[:]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the time delays, the pulsar is now dedispersed for every guessed DM. This is done by using numpy.roll. The threshold for the roll is set by finding the time that each element in the data corresponds to, and moving elements so that they line up with their corresponding elements(come back to give better explanation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dedispersion(DMss, frequency_channels, dispersed_pulse, time_change, n_time, integration_size): \n",
    "    \"\"\"The function takes a set of DMs and dedisperses the given dispersed pulse for all the DMs given.\"\"\"\n",
    "    dedispers = np.zeros((len(DMss),dispersed_pulse.shape[0], dispersed_pulse.shape[1]))  #The dedispersed pulse will have the same shape as the dispersed one\n",
    "    for i in range(len(DMss)):\n",
    "        dedis = dispersed_pulse*0\n",
    "        for j in range(frequency_channels):\n",
    "            dedis[j,:] = np.roll(dispersed_pulse[j,:],-int(round(time_change[i,j]/(time_length/n_time))) )#np.roll shifts \n",
    "                #the data over after the last position and reintroduces it to the first position\n",
    "        dedispers[i,:,:] = dedis\n",
    "    return dedispers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "4.1481e3/1.005"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the above definition, the pulse is dedispersed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dedisperse = dedispersion(DMs, step, abs(dispersed), time_change, n_time, integration_size)  #The dedispersed pulse will have the same shape as the dispersed one\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One DM value for the pulse is then plotted below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "figure()\n",
    "imshow(abs(dedisperse[100]), extent=[np.amin(np.linspace(0,.15, dis_length)),\\\n",
    "                            np.amax(np.linspace(0,.15, dis_length)),\\\n",
    "                            np.amin(pulse_freqs),\\\n",
    "                            np.amax(pulse_freqs)],\\\n",
    "                              aspect='auto')\n",
    "axis_labels('time(s)','Frequency (MHz)',  'Dedispersed Pulsar (DM=600)')\n",
    "colorbar(label = 'power') #This shows the effec the dedispersion had on the 200th DM\n",
    "#savefig(\"generated/Dedispersed_pulsar_DM=600.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To create a DM vs time plot, the array is summed along the frequency axis for every DM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DM_measure = np.sum(dedisperse, axis = 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This new array is then plotted below. Note that the graph appears to have a wave in the region with the highest amplitudes. This is caused by the moire effect which can be shown as the cause with the graph following the one below(see pic for no noise pulsar)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "figure()\n",
    "imshow((abs(DM_measure)**17), extent=[np.amin(np.linspace(0,.15, dis_length)),\\\n",
    "                            np.amax(np.linspace(0,.15, dis_length)),\\\n",
    "                            np.amax(DMs),np.amin(DMs)\n",
    "                            ],\\\n",
    "                              aspect='auto')\n",
    "axis_labels('time(s)','Dispersion Measure',  'Dispersion Measure vs Time of Pulsar (DM=600)')\n",
    "colorbar(label = 'power') #This is the plot of DM vs time\n",
    "#savefig(\"generated/dispersion_measure_vs_time_DM=600_no_noise.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "figure()\n",
    "plot(np.linspace(0,.15, n_time),abs(DM_measure[10])-np.mean(DM_measure[297]), marker='.')\n",
    "axis_labels('time(s)','Amplitude','Amplitude of DM=600')\n",
    "#savefig(\"generated/Amplitude_vs_time_DM=600_no_noise.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the vaience of the non pulse signal needs to be found in order to calculate the signal to noise ratio. This is done by setting a threshold above a certain value, and finding peaks above that threshold. Those peaks are then assumed to belong to signal, and they are removed from the data, in order to find the base noise. The function for this is defined below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pulsefinder(DM_data):\n",
    "    allx_range = []\n",
    "    ally_range = []\n",
    "    for k in range(len(DM_data)):\n",
    "        intensity = DM_data[k] #This defines the input data that we wish to find the peaks in\n",
    "        threshold = np.mean(DM_data[k])+np.std(DM_data[k])#This creates a limit to which peaks with values\n",
    "        #below this number will not be considered in the peak finder\n",
    "        peaksx = [] #This is the values that the peaks are centred around\n",
    "        peaksy = [] #This is the amplitude value of the peaks\n",
    "        peak_dex = (list(intensity)) #This turns the input from an array to a list.\n",
    "        inversep_dex=intensity[::-1] #This reverses the list of the input data\n",
    "        for i in range(len(intensity)-2): #This is the loop that identifies the peaks. two previous values before and after the considered \n",
    "            #point must have values less than the point for said point to be considered a maxima\n",
    "            if intensity[i] >  threshold and intensity[i-2]<intensity[i] and intensity[i-1] <intensity[i] and intensity[i+1] < intensity[i] and intensity[i+2] < intensity[i]:\n",
    "                peaksx.append(np.arange(len(DM_data[k]))[i]) #Here the x value of the peak is appended to a list\n",
    "                peaksy.append(intensity[i]) #The y value of the peak is appended to a list\n",
    "                \n",
    "        #Calculate the FWHM find all point that are part of the peak, or above the median of the data\n",
    "        xrange = []\n",
    "        yrange = []\n",
    "        for i in range(len(peaksy)):\n",
    "            halfmax = peaksy[i]/2 #This caluclates the amplitude of the FWHM\n",
    "            peak_index = peak_dex.index(peaksy[i]) #This identifies the index value of the first peak\n",
    "            pind=(list(inversep_dex)).index(peaksy[i]) #This identifies the index value from the reversed input data array\n",
    "        \n",
    "            for j in range(len((intensity[0:pind]))): #This loop finds the minimum x value that could be considered part of the pulse\n",
    "                if pind == len(intensity)-1: #This ensures that a point that is at the begining is not considered as a peak\n",
    "                    break\n",
    "                ymin = (inversep_dex[pind:][j] - halfmax) #This subtracts the halfmax from the value. Below it is chekced if this difference is below 0\n",
    "                xmins = peak_index-j #This finds the index of the corresponding value\n",
    "                if ymin < 0 or inversep_dex[pind:][j] > peaksy[i] or inversep_dex[pind:][j] <= np.median(intensity) or pind+j+2 > len(intensity) or pind==0: #here the restrictions break the loop if the counter goes below zero, it the yvalue surpasses \n",
    "                 #the yvalue of the peak, or if the yvalue of the index drops below the median of the data set\n",
    "                    break\n",
    "            for j in range(len((intensity[0:peak_index]))):#This loop finds the maximum x value that could be considered part of the pulse\n",
    "                if peak_index == len(intensity)-1: #This ensures that a point that is at the end is not considered as a peak\n",
    "                    break\n",
    "                ymax = (intensity[peak_index:][j] - halfmax) #Does the same as its counterpart ymax above\n",
    "                xmaxs = peak_index+j #Again finds the corresponding index value\n",
    "                if ymax < 0 or intensity[peak_index:][j] > peaksy[i] or intensity[peak_index:][j] < np.median(intensity) or  peak_index+j +2 > len(intensity):\n",
    "                #The above line does the dame as in the previous for loop, and will also break is one of the conditions are met.\n",
    "                    break\n",
    "                \n",
    "        \n",
    "            xrange= xrange+(np.ndarray.tolist(np.arange(len(intensity))[xmins:xmaxs])) #This adds a each new set of values to the existing\n",
    "            #set to make a master set of all values not considered to be noise\n",
    "            yrange= yrange+(np.ndarray.tolist(intensity[xmins:xmaxs])) #This finds the y-value to these points\n",
    "            \n",
    "        allx_range.append(xrange)\n",
    "        ally_range.append(yrange)\n",
    "    return allx_range, ally_range"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function is used below to find the corresponding x and y parameters of the peaks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_ranges,yranges = pulsefinder(DM_measure)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The peaks are then set into a list, and using np.pop, values from the dispersion data are excluded is their x value is included in the list of peaks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_noised_datas = []\n",
    "for k in range(len(DM_measure)):\n",
    "    new_list = np.ndarray.tolist(np.sort(list(set(x_ranges[k]))))\n",
    "    \n",
    "    new_list = new_list+[0]#An extra incedies is added so that when the final useful value is poped below, the loop will\n",
    "    #not throw an error.\n",
    "    noised_data = [] #This is the resulting data points that will be considered noise\n",
    "    for i in np.arange(len(DM_measure[k])):\n",
    "        if i==new_list[0]: #If the counter matches the value of the x_range, the x_range value will be removed, and the \n",
    "            #corresponding data point from the Dispersion data will not be placed in the noise data\n",
    "            new_list.pop(0)\n",
    "            continue\n",
    "        noised_data.append(DM_measure[k][i]) #The values of the DM data that will be considered noise are now added to the noise array\n",
    "    all_noised_datas.append(noised_data)\n",
    "all_noised_datas = np.asarray(all_noised_datas)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "1/2.41"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "figure()\n",
    "plot(all_noised_datas[300], marker='.')\n",
    "plot(all_noised_datas[285], marker='.')\n",
    "\n",
    "#plot(DM_measure[0])\n",
    "axis_labels('time(s)(samples have been cut out, so point locations not accurate)','Amplitude','Noise with a DM of 600, at DM=600')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the new base noise data, the std of each DM data set is found below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now find the standard deviation and the mean of the noise for each DM. More accurate DMs will have a more accurate noise reading.\n",
    "DM_std = np.zeros(len(all_noised_datas))\n",
    "DM_mean = np.zeros(len(all_noised_datas))\n",
    "##Find the mean and std of the noise data\n",
    "for k in range(len(all_noised_datas)):\n",
    "    DM_std[k] = np.std(all_noised_datas[k])\n",
    "    DM_mean[k] = np.mean(all_noised_datas[k])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally to calculate the signal to noise ratio, the peak needs to be calculated. This is done by using those same peaks and averaging them together."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def peakaverage(DM_data,adder):\n",
    "    allx_range = []\n",
    "    ally_range = []\n",
    "    for k in range(len(DM_data)):\n",
    "        intensity = DM_data[k] #This defines the input data that we wish to find the peaks in\n",
    "        threshold = np.mean(DM_data[k])+np.std(DM_data[k])#This creates a limit to which peaks with values\n",
    "        #below this number will not be considered in the peak finder\n",
    "        peaksx = [] #This is the values that the peaks are centred around\n",
    "        peaksy = [] #This is the amplitude value of the peaks\n",
    "        peak_dex = (list(intensity)) #This turns the input from an array to a list.\n",
    "        inversep_dex=intensity[::-1] #This reverses the list of the input data\n",
    "        for i in range(len(intensity)-2): #This is the loop that identifies the peaks. two previous values before and after the considered \n",
    "            #point must have values less than the point for said point to be considered a maxima\n",
    "            if intensity[i] >  threshold and intensity[i-2]<intensity[i] and intensity[i-1] <intensity[i] and intensity[i+1] < intensity[i] and intensity[i+2] < intensity[i]:\n",
    "                peaksx.append(np.arange(len(DM_data[k]))[i]) #Here the x value of the peak is appended to a list\n",
    "                peaksy.append(intensity[i]) #The y value of the peak is appended to a list\n",
    "        xrange = []\n",
    "        yrange = []\n",
    "        for i in range(len(peaksy)):\n",
    "            halfmax = peaksy[i]/2 #This caluclates the amplitude of the FWHM\n",
    "            peak_index = peak_dex.index(peaksy[i]) #This identifies the index value of the first peak\n",
    "            pind=(list(inversep_dex)).index(peaksy[i]) #This identifies the index value from the reversed input data array\n",
    "        \n",
    "            for j in range(len((intensity[0:pind]))): #This loop finds the minimum x value that could be considered part of the pulse\n",
    "                if pind == len(intensity)-1: #This ensures that a point that is at the begining is not considered as a peak\n",
    "                    break\n",
    "                ymin = (inversep_dex[pind:][j] - halfmax) #This subtracts the halfmax from the value. Below it is chekced if this difference is below 0\n",
    "                xmins = peak_index-j #This finds the index of the corresponding value\n",
    "                if ymin < 0 or inversep_dex[pind:][j] > peaksy[i] or pind+j+2 > len(intensity) or pind==0: #here the restrictions break the loop if the counter goes below zero, it the yvalue surpasses \n",
    "                 #the yvalue of the peak, or if the yvalue of the index drops below the median of the data set\n",
    "                    break\n",
    "            for j in range(len((intensity[0:peak_index]))):#This loop finds the maximum x value that could be considered part of the pulse\n",
    "                if peak_index == len(intensity)-1: #This ensures that a point that is at the end is not considered as a peak\n",
    "                    break\n",
    "                ymax = (intensity[peak_index:][j] - halfmax) #Does the same as its counterpart ymax above\n",
    "                xmaxs = peak_index+j #Again finds the corresponding index value\n",
    "                if ymax < 0 or intensity[peak_index:][j] > peaksy[i] or  peak_index+j +2 > len(intensity):\n",
    "                #The above line does the dame as in the previous for loop, and will also break is one of the conditions are met.\n",
    "                    break\n",
    "                \n",
    "        \n",
    "            xrange= xrange+(np.ndarray.tolist(np.arange(len(intensity))[xmins:xmaxs])) #This adds a each new set of values to the existing\n",
    "            #set to make a master set of all values not considered to be noise\n",
    "            yrange= yrange+(np.ndarray.tolist(intensity[xmins:xmaxs])) #This finds the y-value to these points\n",
    "            \n",
    "        allx_range.append(xrange)\n",
    "        ally_range.append(yrange)\n",
    "    \n",
    "    peaks_av = np.zeros(len(ally_range))\n",
    "    for k in range(len(ally_range)):\n",
    "        peaks_av[k] = np.mean(ally_range[k]+np.mean(adder[k]))\n",
    "    \n",
    "    return peaks_av"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DM_meaned = np.zeros((len(DM_measure),n_time))\n",
    "for i in range(len(DM_measure)):\n",
    "    DM_meaned[i] = DM_measure[i]-np.mean(DM_measure[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ypeaks_av = peakaverage(DM_meaned, DM_measure)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "figure()\n",
    "plot((ypeaks_av[:]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Another method to find appropriate DMs is to use matched filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DM_measure[200].max()-DM_measure[201].max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(DM_mean[295])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(ypeaks_av[249:350]-DM_mean[249:350])/DM_std[249:350]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "1500000/400"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}