40m QIL Cryo_Lab CTN SUS_Lab TCS_Lab OMC_Lab CRIME_Lab FEA ENG_Labs OptContFac Mariner WBEEShop
  40m Log, Page 34 of 327  Not logged in ELOG logo
ID Date Authorup Type Category Subject
  16151   Fri May 21 09:44:52 2021 Ian MacMillanUpdateCDSSUS simPlant model

The transfer function given in the previous post was slightly incorrect the units did not make sense the new function is:

\frac{x}{F}=\frac{1}{m\omega_0^2-m\omega^2+im\frac{\omega_0 \omega }{Q}}

I have attached a quick derivation below in attachment 1

Attachment 1: Transfer_Function_of_Damped_Harmonic_Oscillator.pdf
Transfer_Function_of_Damped_Harmonic_Oscillator.pdf
  16153   Fri May 21 14:36:20 2021 Ian MacMillanUpdateCDSSUS simPlant model

The plant transfer function of the pendulum in the s domain is:

H(s)=\frac{x(s)}{F(s)}=\frac{1}{ms^2+m\frac{\omega_0}{Q}s+m\omega_0^2}

Using Foton to make a plot of the TF needed and using m=40kg, w0=3Hz, and Q=50 (See attachment 1). It is easiest to enter the above filter using RPoly and saved it as Plant_V1

Attachment 1: Plant_Mod_TF.pdf
Plant_Mod_TF.pdf
  16177   Thu Jun 3 13:06:47 2021 Ian MacMillanUpdateCDSSUS simPlant model

I was able to measure the transfer function of the plant filter module from the channel X1:SUP-C1_SUS_SINGLE_PLANT_Plant_POS_Mod_EXC to X1:SUP-C1_SUS_SINGLE_PLANT_Plant_POS_Mod_OUT. The resulting transfer function is shown below. I have also attached the raw data for making the graph.

Next, I will make a script that will make the photon filters for all the degrees of freedom and start working on the matrix version of the filter module so that there can be multiple degrees of freedom.

Attachment 1: SingleSusPlantTF.pdf
SingleSusPlantTF.pdf
Attachment 2: SUS_PLANT_TF.zip
  16191   Mon Jun 7 17:49:19 2021 Ian MacMillanUpdateCDSSUS simPlant model

Added difference to the graph. I included the code so that others could see what it looks like and use it for easy use.

Attachment 1: SingleSusPlantTF.pdf
SingleSusPlantTF.pdf
Attachment 2: TF_Graph_Code.zip
  16195   Wed Jun 9 13:50:48 2021 Ian MacMillanUpdateCDSSUS simPlant model

I have attached an updated transfer function graph with the residual easier to see. I thought here I would include a better explanation of what this transfer function was measuring.

This transfer function was mainly about learning how to use DTT and Foton to make and measure transfer functions. Therefore it is just measuring across a single CDS filter block. X1SUP_ISOLATED_C1_SUS_SINGLE_PLANT_Plant_POS_Mod block to be specific. This measurement shows that the block is doing what I programmed it to do with Foton. The residual is probably just because the measured TF had fewer points than the calculated one.

The next step is to take a closed-loop TF of the system and the control module.

After that, I want to add more degrees of freedom to the model. both in the plant and in the controls.

Attachment 1: SingleSusPlantTF.pdf
SingleSusPlantTF.pdf
  16201   Tue Jun 15 11:46:40 2021 Ian MacMillanUpdateCDSSUS simPlant model

I have added more degrees of freedom. The model includes x, y, z, pitch, yaw, roll and is controlled by a matrix of transfer functions (See Attachment 2). I have added 5 control filters to individually control UL, UR, LL, LR, and side. Eventually, this should become a matrix too but for the moment this is fine.

Note the Unit delay blocks in the control in Attachment 1. The model will not compile without these blocks.

Attachment 1: x1sup_isolated-6-15-v1.pdf
x1sup_isolated-6-15-v1.pdf
Attachment 2: C1_SUS_SINGLE_PLANT-6-15-v1.pdf
C1_SUS_SINGLE_PLANT-6-15-v1.pdf
  16208   Thu Jun 17 11:19:37 2021 Ian MacMillanUpdateCDSCDS Upgrade

Jon and I tested the ADC and DAC cards in both of the systems on the test stand. We had to swap out an 18-bit DAC for a 16-bit one that worked but now both machines have at least one working ADC and DAC.

[Still working on this post. I need to look at what is in the machines to say everything ]

  16217   Mon Jun 21 17:15:49 2021 Ian MacMillanUpdateCDSCDS Upgrade

Anchal and I wrote a script (Attachment 1) that will test the ADC and DAC connections with inputs on the INMON from -3000 to 3000. We could not run it because some of the channels seemed to be frozen. 

Attachment 1: DAC2ADC_Test.py
´╗┐import os
import time
import numpy as np
import subprocess
from traceback import print_exc
import argparse


def grabInputArgs():
    parser = argparse.ArgumentParser(
... 75 more lines ...
  16220   Tue Jun 22 16:53:01 2021 Ian MacMillanUpdateCDSFront-End Assembly and Testing

The channels on both the C1BHD and C1SUS2 seem to be frozen: they arent updating and are holding one value. To fix this Anchal and I tried:

  • restarting the computers 
    • restarting basically everything including the models
  • Changing the matrix values
  • adding filters
  • messing with the offset 
  • restarting the network ports (Paco suggested this apparently it worked for him at some point)
  • Checking to make sure everything was still connected inside the case (DAC, ADC, etc..)

I wonder if Jon has any ideas. 

  16224   Thu Jun 24 17:32:52 2021 Ian MacMillanUpdateCDSFront-End Assembly and Testing

Anchal and I ran tests on the two systems (C1-SUS2 and C1-BHD). Attached are the results and the code and data to recreate them.

We connected one DAC channel to one ADC channel and thus all of the results represent a DAC/ADC pair. We then set the offset to different values from -3000 to 3000 and recorded the measured signal. I then plotted the response curve of every DAC/ADC pair so each was tested at least once.

There are two types of plots included in the attachments

1) a summary plot found on the last pages of the pdf files. This is a quick and dirty way to see if all of the channels are working. It is NOT a replacement for the other plots. It shows all the data quickly but sacrifices precision.

2) In an in-depth look at an ADC/DAC pair. Here I show the measured value for a defined DC offset. The Gain of the system should be 0.5 (put in an offset of 100 and measure 50). I included a line to show where this should be. I also plotted the difference between the 0.5 gain line and the measured data. 

As seen in the provided plots the channels get saturated after about the -2000 to 2000 mark, which is why the difference graph is only concentrated on -2000 to 2000 range. 

Summary: all the channels look to be working they all report very little deviation off of the theoretical gain. 

Note: ADC channel 31 is the timing signal so it is the only channel that is wildly off. It is not a measurement channel and we just measured it by mistake.

Attachment 1: C1-SU2_Channel_Responses.pdf
C1-SU2_Channel_Responses.pdf C1-SU2_Channel_Responses.pdf C1-SU2_Channel_Responses.pdf C1-SU2_Channel_Responses.pdf C1-SU2_Channel_Responses.pdf C1-SU2_Channel_Responses.pdf C1-SU2_Channel_Responses.pdf C1-SU2_Channel_Responses.pdf
Attachment 2: C1-BHD_Channel_Responses.pdf
C1-BHD_Channel_Responses.pdf C1-BHD_Channel_Responses.pdf C1-BHD_Channel_Responses.pdf C1-BHD_Channel_Responses.pdf C1-BHD_Channel_Responses.pdf C1-BHD_Channel_Responses.pdf C1-BHD_Channel_Responses.pdf C1-BHD_Channel_Responses.pdf
Attachment 3: CDS_Channel_Test.zip
  16230   Wed Jun 30 14:09:26 2021 Ian MacMillanUpdateCDSSUS simPlant model

I have looked at my code from the previous plot of the transfer function and realized that there is a slight error that must be fixed before we can analyze the difference between the theoretical transfer function and the measured transfer function.

The theoretical transfer function, which was generated from Photon has approximately 1000 data points while the measured one has about 120. There are no points between the two datasets that have the same frequency values, so they are not directly comparable. In order to compare them I must infer the data between the points. In the previous post [16195] I expanded the measured dataset. In other words: I filled in the space between points linearly so that I could compare the two data sets. Using this code:

#make values for the comparison
tck_mag = splrep(tst_f, tst_mag) # get bspline representation given (x,y) values
gen_mag = splev(sim_f, tck_mag) # generate intermediate values
dif_mag=[]
for x in range(len(gen_mag)):
    dif_mag.append(gen_mag[x]-sim_mag[x]) # measured minus predicted

tck_ph = splrep(tst_f, tst_ph) # get bspline representation given (x,y) values
gen_ph = splev(sim_f, tck_ph) # generate intermediate values
dif_ph=[]
for x in range(len(gen_ph)):
    dif_ph.append(gen_ph[x]-sim_ph[x])

At points like a sharp peak where the measured data set was sparse compared to the peak, the difference would see the difference between the intermediate “measured” values and the theoretical ones, which would make the difference much higher than it really was.

To fix this I changed the code to generate the intermediate values for the theoretical data set. Using the code here:

tck_mag = splrep(sim_f, sim_mag) # get bspline representation given (x,y) values
gen_mag = splev(tst_f, tck_mag) # generate intermediate values
dif_mag=[]
for x in range(len(tst_mag)):
    dif_mag.append(tst_mag[x]-gen_mag[x])#measured minus predicted

tck_ph = splrep(sim_f, sim_ph) # get bspline representation given (x,y) values
gen_ph = splev(tst_f, tck_ph) # generate intermediate values
dif_ph=[]
for x in range(len(tst_ph)):
    dif_ph.append(tst_ph[x]-gen_ph[x])

Because this dataset has far more values (about 10 times more) the previous problem is not such an issue. In addition, there is never an inferred measured value used. That makes it more representative of the true accuracy of the real transfer function.

This is an update to a previous plot, so I am still using the same data just changing the way it is coded. This plot/data does not have a Q of 1000. That plot will be in a later post along with the error estimation that we talked about in this week’s meeting.

The new plot is shown below in attachment 1. Data and code are contained in attachment 2

Attachment 1: SingleSusPlantTF.pdf
SingleSusPlantTF.pdf
Attachment 2: Plant_TF_Test.zip
  16289   Mon Aug 23 15:25:59 2021 Ian MacMillanUpdateCDSSUS simPlant model

I am adding a State-space block to the SimPlant cds model using the example Chris gave. I made a new folder in controls called SimPlantStateSpace. wI used the code below to make a state-space LTI model with a 1D pendulum then I converted it to a discrete system using c2d matlab function. Then I used these in the rtss.m file to create the state space code I need in the SimPlantStateSpace_1D_model.h file.

%sys_model.m

Q = 1000;
phi = 1/Q;
g = 9.806;
m = 0.24; % mass of pendulum
l = 0.248; %length of pendulum
w_0 = sqrt(g/l);

f=16000 %this is the frequency of the channel that will be used

A = [0 1; -w_0^2*(1+1/Q*1i) -w_0/Q]
B = [0; 1/m];
C = [1 0];
D = [0];
sys_dc = ss(A,B,C,D)

sys=c2d(sys_dc, 1/f)

This code outputs the discrete state space that is added to the header file attached.

Attachment 1: SimPlantStateSpace.zip
  16339   Thu Sep 16 14:08:14 2021 Ian MacMillanFrogs Tour

I gave some of the data analysts a look around because they asked and nothing was currently going on in the 40m. Nothing was changed.

  16348   Mon Sep 20 15:42:44 2021 Ian MacMillanSummaryComputersQuantization Code Summary

This post serves as a summary and description of code to run to test the impacts of quantization noise on a state-space implementation of the suspension model.

Purpose: We want to use a state-space model in our suspension plant code. Before we can do this we want to test to see if the state-space model is prone to problems with quantization noise. We will compare two models one for a standard direct-ii filter and one with a state-space model and then compare the noise from both. 

Signal Generation:

First I built a basic signal generator that can produce a sine wave for a specified amount of time then can produce a zero signal for a specified amount of time. This will let the model ring up with the sine wave then decay away with the zero signal. This input signal is generated at a sample rate of 2^16 samples per second then stored in a numpy array. I later feed this back into both models and record their results.

State-space Model:

The code can be seen here

The state-space model takes in the list of excitation values and feeds them through a loop that calculates the next value in the output.

Given that the state-space model follows the form

  \dot{x}(t)=\textbf{A}x(t)+ \textbf{B}u(t)   and  y(t)=\textbf{C}x(t)+ \textbf{D}u(t) ,

the model has three parts the first equation, an integration, and the second equation.

  1. The first equation takes the input x and the excitation u and generates the x dot vector shown on the left-hand side of the first state-space equation.
  2. The second part must integrate x to obtain the x that is found in the next equation. This uses the velocity and acceleration to integrate to the next x that will be plugged into the second equation
  3. The second equation in the state space representation takes the x vector we just calculated and then multiplies it with the sensing matrix C. we don't have a D matrix so this gives us the next output in our system

This system is the coded form of the block diagram of the state space representation shown in attachment 1

Direct-II Model:

The direct form 2 filter works in a much simpler way. because it involves no integration and follows the block diagram shown in Attachment 2, we can use a single difference equation to find the next output. However, the only complication that comes into play is that we also have to keep track of the w(n) seen in the middle of the block diagram. We use these two equations to calculate the output value

y[n]=b_0 \omega [n]+b_1 \omega [n-1] +b_2 \omega [n-2],  where w[n] is  \omega[n]=x[n] - a_1 \omega [n-1] -a_2 \omega[n-2]

Bit length Control:

To control the bit length of each of the models I typecast all the inputs using np.float then the bit length that I want. This simulates the computer using only a specified bit length. I have to go through the code and force the code to use 128 bit by default. Currently, the default is 64 bit which so at the moment I am limited to 64 bit for highest bit length. I also need to go in and examine how numpy truncates floats to make sure it isn't doing anything unexpected.

Bode Plot: 

The bode plot at the bottom shows the transfer function for both the IIR model and the state-space model. I generated about 100 seconds of white noise then computed the transfer function as 

G(f) = \frac{P_{csd}(f)}{P_{psd}(f)}

which is the cross-spectral density divided by the power spectral density. We can see that they match pretty closely at 64 bits. The IIR direct II model seems to have more noise on the surface but we are going to examine that in the next elog

 

Attachment 1: 472px-Typical_State_Space_model.svg.png
472px-Typical_State_Space_model.svg.png
Attachment 2: Biquad_filter_DF-IIx.svg.png
Biquad_filter_DF-IIx.svg.png
Attachment 3: SS-IIR-TF.pdf
SS-IIR-TF.pdf
  16355   Wed Sep 22 14:22:35 2021 Ian MacMillanSummaryComputersQuantization Noise Calculation Summary

Now that we have a model of how the SS and IIR filters work we can get to the problem of how to measure the quantization noise in each of the systems. Den Martynov's thesis talks a little about this. from my understanding: He measured quantization noise by having two filters using two types of variables with different numbers of bits. He had one filter with many more bits than the second one. He fed the same input signal to both filters then recorded their outputs x_1 and x_2, where x_2 had the higher number of bits. He then took the difference x_1-x_2. Since the CDS system uses double format, he assumes that quantization noise scales with mantissa length. He can therefore extrapolate the quantization noise for any mantissa length.

Here is the Code that follows the following procedure (as of today at least)

This problem is a little harder than I had originally thought. I took Rana's advice and asked Aaron about how he had tackled a similar problem. We came up with a procedure explained below (though any mistakes are my own):

  1. Feed different white noise data into three of the same filter this should yield the following equation: \textbf{S}_i^2 =\textbf{S}_{ni}^2+\textbf{S}_x^2, where \textbf{S}_i^2 is the power spectrum of the output for the ith filter,  \textbf{S}_{ni}^2 is the noise filtered through an "ideal" filter with no quantization noise, and  \textbf{S}_x^2 is the power spectrum of the quantization noise. Since we are feeding random noise into the input the power of the quantization noise should be the same for all three of our runs.
  2. Next, we have our three outputs:  \textbf{S}_1^2,  \textbf{S}_2^2, and  \textbf{S}_3^2 that follow the equations: 

\textbf{S}_1^2 =\textbf{S}_{n1}^2+\textbf{S}_x^2

\textbf{S}_2^2 =\textbf{S}_{n2}^2+\textbf{S}_x^2

\textbf{S}_3^2 =\textbf{S}_{n3}^2+\textbf{S}_x^2

From these three equations, we calculate the three quantities: \textbf{S}_{12}^2\textbf{S}_{23}^2, and \textbf{S}_{13}^2 which are calculated by:

\textbf{S}_{12}^2 =\textbf{S}_{1}^2-\textbf{S}_2^2\approx \textbf{S}_{n1}^2 -\textbf{S}_{n2}^2

\textbf{S}_{23}^2 =\textbf{S}_{2}^2-\textbf{S}_3^2\approx \textbf{S}_{n2}^2 -\textbf{S}_{n3}^2

\textbf{S}_{13}^2 =\textbf{S}_{1}^2-\textbf{S}_3^2\approx \textbf{S}_{n1}^2 -\textbf{S}_{n3}^2

from these quantities, we can calculate three values: \bar{\textbf{S}}_{n1}^2\bar{\textbf{S}}_{n2}^2, and \bar{\textbf{S}}_{n3}^2 since these are just estimates we are using a bar on top. These are calculated using:

\bar{\textbf{S}}_{n1}^2\approx\frac{1}{2}(\textbf{S}_{12}^2+\textbf{S}_{13}^2+\textbf{S}_{23}^2)

\bar{\textbf{S}}_{n2}^2\approx\frac{1}{2}(-\textbf{S}_{12}^2+\textbf{S}_{13}^2+\textbf{S}_{23}^2)

\bar{\textbf{S}}_{n3}^2\approx\frac{1}{2}(\textbf{S}_{12}^2+\textbf{S}_{13}^2-\textbf{S}_{23}^2)

using these estimates we can then estimate  \textbf{S}_{x}^2  using the formula:

\textbf{S}_{x}^2 \approx \textbf{S}_{1}^2 - \bar{\textbf{S}}_{n1}^2 \approx \textbf{S}_{2}^2 - \bar{\textbf{S}}_{2}^2 \approx \textbf{S}_{3}^2 - \bar{\textbf{S}}_{n3}^2

we can average the three estimates for  \textbf{S}_{x}^2  to come up with one estimate.

This procedure should be able to give us a good estimate of the quantization noise. However, in the graph shown in the attachments below show that the noise follows the transfer function of the model to begin with. I would not expect this to be true so I believe that there is an error in the above procedure or in my code that I am working on finding. I may have to rework this three-corner hat approach. I may have a mistake in my code that I will have to go through.

I would expect the quantization noise to be flatter and not follow the shape of the transfer function of the model. Instead, we have what looks like just the result of random noise being filtered through the model. 

Next steps:

The first real step is being able to quantify the quantization noise but after I fix the issues in my code I will be able to start liking at optimal model design for both the state-space model and the direct form II model. I have been looking through the book "Quantization noise" by Bernard Widrow and Istvan Kollar which offers some good insights on how to minimize quantization noise. 

Attachment 1: IIR64-bitnoisespectrum.pdf
IIR64-bitnoisespectrum.pdf
  16360   Mon Sep 27 12:12:15 2021 Ian MacMillanSummaryComputersQuantization Noise Calculation Summary

I have not been able to figure out a way to make the system that Aaron and I talked about. I'm not even sure it is possible to pull the information out of the information I have in this way. Even the book uses a comparison to a high precision filter as a way to calculate the quantization noise:

"Quantization noise in digital filters can be studied in simulation by comparing the behavior of the actual quantized digital filter with that of a refrence digital filter having the same structure but whose numerical calculations are done extremely accurately."
-Quantization Noise by Bernard Widrow and Istvan Kollar (pg. 416)

Thus I will use a technique closer to that used in Den Martynov's thesis (see appendix B starting on page 171). A summary of my understanding of his method is given here:

A filter is given raw unfiltered gaussian data f(t) then it is filtered and the result is the filtered data x(t) thus we get the result: f(t)\rightarrow x(t)=x_N(t)+x_q(t)  where x_N(t) is the raw noise filtered through an ideal filter and x_q(t) is the difference which in this case is the quantization noise. Thus I will input about 100-1000 seconds of the same white noise into a 32-bit and a 64-bit filter. (hopefully, I can increase the more precise one to 128 bit in the future) then I record their outputs and subtract the from each other. this should give us the Quantization error e(t):
e(t)=x_{32}(t)-x_{64}(t)=x_{N_{32}}(t)+x_{q_{32}}(t) - x_{N_{64}}(t)-x_{q_{64}}(t)
and since x_{N_{32}}(t)=x_{N_{64}}(t) because they are both running through ideal filters:
e(t)=x_{N}(t)+x_{q_{32}}(t) - x_{N}(t)-x_{q_{64}}(t)
e(t)=x_{q_{32}}(t) -x_{q_{64}}(t)
and since in this case, we are assuming that the higher bit-rate process is essentially noiseless we get the Quantization noise x_{q_{32}}(t).

If we make some assumptions, then we can actually calculate a more precise version of the quantization noise:

"Since aLIGO CDS system uses double precision format, quantization noise is extrapolated assuming that it scales with mantissa length"
-Denis Martynov's Thesis (pg. 173)

From this assumption, we can say that the noise difference between the 32-bit and 64-bit filter outputs:  x_{q_{32}}(t)-x_{q_{64}}(t)  is proportional to the difference between their mantissa length. by averaging over many different bit lengths, we can estimate a better quantization noise number.

I am building the code to do this in this file

  16361   Mon Sep 27 16:03:15 2021 Ian MacMillanSummaryComputersQuantization Noise Calculation Summary

I have coded up the procedure in the previous post: The result does not look like what I would expect. 

As shown in Attachment1 I have the power spectrum of the 32-bit output and the 64-bit output as well as the power spectrum of the two subtracted time series as well as the subtracted power spectra of both. unfortunately, all of them follow the same general shape of the raw output of the filter. 

I would not expect quantization noise to follow the shape of the filter. I would instead expect it to be more uniform. If anything I would expect the quantization noise to increase with frequency. If a high-frequency signal is run through a filter that has high quantization noise then it will severely degrade: i.e. higher quantization noise. 

This is one reason why I am confused by what I am seeing here. In all cases including feeding the same and different white noise into both filters, I have found that the calculated quantization noise is proportional to the response of the filter. this seems wrong to me so I will continue to play around with it to see if I can gain any intuition about what might be happening.

Attachment 1: DeltaNoiseSpectrum.pdf
DeltaNoiseSpectrum.pdf
  16366   Thu Sep 30 11:46:33 2021 Ian MacMillanSummaryComputersQuantization Noise Calculation Summary

First and foremost I have the updated bode plot with the mode moved to 10 Hz. See Attachment 1. Note that the comparison measurement is a % difference whereas in the previous bode plot it was just the difference. I also wrapped the phase so that jumps from -180 to 180 are moved down. This eliminates massive jumps in the % difference. 

Next, I have two comparison plots: 32 bit and 64bit. As mentioned above I moved the mode to 10 Hz and just excited both systems at 3.4283Hz with an amplitude of 1. As we can see on the plot the two models are practically the same when using 64bits. With the 32bit system, we can see that the noise in the IIR filter is much greater than in the State-space model at frequencies greater than our mode.

Note about windowing and averaging: I used a Hanning window with averaging over 4 neighbor points. I came to this number after looking at the results with less averaging and more averaging. In the code, this can be seen as nperseg=num_samples/4 which is then fed into signal.welch

Attachment 1: SS-IIR-Bode.pdf
SS-IIR-Bode.pdf
Attachment 2: PSD_32bit.pdf
PSD_32bit.pdf
Attachment 3: PSD_64bit.pdf
PSD_64bit.pdf
  16403   Thu Oct 14 16:38:26 2021 Ian MacMillanUpdateGeneralKicking optics in freeSwing measurment

[Ian, Anchal]

We are going to kick the optics tonight at 2am.

The optics we will kick are the PRM BS ITMX ITMY ETMX ETMY

We will kick each one once and record for 2000 seconds and the log files will be placed in users/ian/20211015_FreeSwingTest/logs.

  16406   Fri Oct 15 12:14:27 2021 Ian MacMillanUpdateGeneralKicking optics in freeSwing measurment

[Ian, Anchal]

we ran the free swinging test last night and the results match up with in 1/10th of a Hz. We calculated the peak using the getPeakFreqs2 script to find the peaks and they are close to previous values from 2016.

In attachment 1 you will see the results of the test for each optic.

The peak values are as follows:

Optic POS (Hz) PIT (Hz) YAW (Hz) SIDE (Hz)
PRM 0.94 0.96 0.99 0.99
MC2 0.97 0.75 0.82 0.99
ETMY 0.98 0.98 0.95 0.95
MC1 0.97 0.68 0.80 1.00
ITMX 0.95 0.68 0.68 0.98
ETMX 0.96 0.73 0.85 1.00
BS 0.99 0.74 0.80 0.96
ITMY 0.98 0.72 0.72 0.98
MC3 0.98 0.77 0.84 0.97

The results from 2016 can be found at: /cvs/cds/rtcdt/caltech/c1/scripts/SUS/PeakFit/parameters2.m

Attachment 1: 20211015_Kicktest_plot.pdf
20211015_Kicktest_plot.pdf 20211015_Kicktest_plot.pdf 20211015_Kicktest_plot.pdf 20211015_Kicktest_plot.pdf 20211015_Kicktest_plot.pdf 20211015_Kicktest_plot.pdf 20211015_Kicktest_plot.pdf 20211015_Kicktest_plot.pdf
  16414   Tue Oct 19 18:20:33 2021 Ian MacMillanSummaryCDSc1sus2 DAC to ADC test

I ran a DAC to ADC test on c1sus2 channels where I hooked up the outputs on the DAC to the input channels on the ADC. We used different combinations of ADCs and DACs to make sure that there were no errors that cancel each other out in the end. I took a transfer function across these channel combinations to reproduce figure 1 in T2000188.

As seen in the two attached PDFs the channels seem to be working properly they have a flat response with a gain of 0.5 (-6 dB). This is the response that is expected and is the result of the DAC signal being sent as a single ended signal and the ADC receiving as a differential input signal. This should result in a recorded signal of 0.5 the amplitude of the actual output signal.

The drop off on the high frequency end is the result of the anti-aliasing filter and the anti-imaging filter. Both of these are 8-pole elliptical filters so when combined we should get a drop off of 320dB per decade. I measured the slope on the last few points of each filter and the averaged value was around 347dB per decade. This is slightly steeper than expected but since it is to cut off higher frequencies it shouldn't have an effect on the operation of the system. Also it is very close to the expected value.

The ripples seen before the drop off are also an effect of the elliptical filters and are seen in T2000188.

Note: the transfer function that doesn't seem to match the others is the heartbeat timing signal.

Attachment 1: data3_Plots.pdf
data3_Plots.pdf data3_Plots.pdf data3_Plots.pdf data3_Plots.pdf data3_Plots.pdf data3_Plots.pdf data3_Plots.pdf data3_Plots.pdf
Attachment 2: data2_Plots.pdf
data2_Plots.pdf data2_Plots.pdf data2_Plots.pdf data2_Plots.pdf data2_Plots.pdf data2_Plots.pdf data2_Plots.pdf data2_Plots.pdf
  16423   Fri Oct 22 17:35:08 2021 Ian MacMillanSummaryPEMParticle counter setup near BS Chamber

I have done some reading about where would be the best place to put the particle counter. The ISO standard (14644-1:2015) for cleanrooms is one every 1000 m^2 so one for every 30m x 30m space. We should have the particle counter reasonably close to the open chamber and all the manufactures that I read about suggest a little more than 1 every 30x30m. We will have it much closer than this so it is nice to know that it should still get a good reading. They also suggest keeping it in the open and not tucked away which is a little obvious. I think the best spot is attached to the cable tray that is right above the door to the control room. This should put it out of the way and within about 5m of where we are working. I ordered some cables to route it over there last night so when they come in I can put it up there. 

  16430   Tue Oct 26 18:24:00 2021 Ian MacMillanSummaryCDSc1sus2 DAC to ADC test

[Ian, Anchal, Paco]

After the Koji found that there was a problem with the power source Anchal and I fixed the power then reran the measurment. The only change this time around is that I increased the excitation amplitude to 100. In the first run the excitation amplitude was 1 which seemed to come out noise free but is too low to give a reliable value.

link to previous results

The new plots are attached.

Attachment 1: data2_Plots.pdf
data2_Plots.pdf data2_Plots.pdf data2_Plots.pdf data2_Plots.pdf data2_Plots.pdf data2_Plots.pdf data2_Plots.pdf data2_Plots.pdf
Attachment 2: data3_Plots.pdf
data3_Plots.pdf data3_Plots.pdf data3_Plots.pdf data3_Plots.pdf data3_Plots.pdf data3_Plots.pdf data3_Plots.pdf data3_Plots.pdf
  16447   Wed Nov 3 16:55:23 2021 Ian MacMillanSummarySUSSUS Plant Plan for New Optics

[Ian, Tega, Raj]

This is the rough plan for the testing of the new suspension models with the created plant model. We will test the suspensions on the plant model before we implement them into the full

  • Get State-space matrices from the surf model for the SOS. Set up simplant model on teststand
    • The state-space model is only 3 degrees of freedom. (even the surf's model)
    • There are filter modules that have the 6 degrees of freedom for the suspensions. We will use these instead. I have implemented them in the same suspension model that would hold the state space model. If we ever get the state space matrices then we can easily substitute them.
  • Load new controller model onto test stand. This new controller will be a copy of an existing suspension controller.
  • Hook up controller to simplant.  These should form a closed loop where the outputs from the controller go into the plant inputs and the plant outputs go to the controller inputs.
  • Do tests on set up.
    • Look at the step response for each degree of freedom. Then compare them to the results from an existing optic. 
      • Also, working with Raj let him do the same model in python then compare the two.
  • Make sure that the data is being written to the local frame handler.

MEDM file location

/opt/rtcds/userapps/release/sus/common/medm/hsss_tega_gautam

run using 

For ITMX display, use:

hsss_tega_gautam>medm -x -macro "site=caltech,ifo=c1,IFO=C1,OPTIC=ITMX,SUSTYPE=IM,DCU_ID=21,FEC=45" SUS_CUST_HSSS_OVERVIEW.adl

  16457   Mon Nov 8 17:52:22 2021 Ian MacMillanUpdateSUSSetting up suspension test model

[Ian, Tega]

We combined a controler and a plant model into a single modle (See first attachment) called x1sus_cp.mdl in the userapps folder of the cymac in c1sim. This model combines 2 blocks: the controler block which is used to control the current optics and is found in cvs/cds/rtcds/userapps/release/sus/c1/models/c1sus.mdl further the control block we are using comes from the same path but from the c1sup.mdl model. This plant model is the bases for all of my custom plant models and thus is a good starting point for the testing. It is also ideal because I know it can beeasily altered into a my state-space plant model. However, we had to make a few adjustments to get the model up to date for the cds system. So it is now a unique block.

These two library blocks are set in the userapps/lib folder on the cymac. This is the lib file that the docker system looks to when it is compiling models. For a quick overview see this. All other models have been removed from the MatLab path so that when we open x1sus_cp.mdl in MatLab it is using the same models it will compile with.

We could not find the rtbitget library part, but chris pointed us to userapps, and we copied it over using: scp /opt/rtcds/userapps/trunk/cds/common/models/rtbitget.mdl controls@c1sim:/home/controls/simLink/lib.

NOTE TO FUTURE IAN: don't forget that unit delays exist.

Next step: now that we have a model that is compiling and familiar we need to make medm screens. We will use the auto mdl2adl for this so that it is quick. Then we can start adding our custom pieces one by one so that we know that they are working. We will also work with Raj to get an independent python model working. Which will allow us to compare the cds and python models.

Attachment 1: x1sus_cp.png
x1sus_cp.png
  16460   Tue Nov 9 13:40:02 2021 Ian MacMillanSummaryComputer Scripts / ProgramsSUS Plant Plan for New Optics

[Ian, Tega]

After talking with Rana we have an updated plan. We will be working on this plan step by step in this order.

  1. Remove c1sim from the test stand rack and move it to the rack in the office next to the printer. When connecting it we will NOT connect it to the Martian network! This is to make sure that nothing is connected to the 40m system and we can't mess anything up.
  2. Once we have moved the computer over physically, we will need to update anyone who uses it on how to connect to it. The way we connect to it will have changed.
  3. Now that we have the computer moved and everyone can connect to it we will work on the model. Currently, we have the empty models connected.
    1. recompile the model since we moved the computer.
    2. verify that nothing has changed in the move and the model can still operate and compile properly
  4. The model has the proper structure but we need to fill it with the proper filters and such
    1. For the Plant model
      1. To get it up and running quickly we will use the premade plant filters for the plant model. These filters were made for the c1sup.mdl and should work in our modified plant model. This will allow us to verify that everything is working. And allow us to run tests on the system.
      2. We need to update the model and add the state space block. (we are skipping this step for now because we are fast-tracking the testing)  
        1. Check with Chris to make sure that this is the right way to do it. I am pretty sure it is, but I don't know anything
        2. Make the 6 DOF state-space matrix. We only have a three DOF one. The surf never made a 6 DOF. 
        3. Make the block to input into the model
        4. make a switch that will allow us to switch between the state-space model and the filter block
    2. For the controller
      1. Load filter coefficients for the controller model from one of the current optics and use this as a starting point.
      2. Add medm screens for the controller and plant. We are skipping this for now because we want results and we don't care if the screens look nice and are useable at the moment.
  5.  Test the model
    1. we will take an open-loop transfer function of all six of the DOFs to all other DOFs which will leave us with 36 TFs. Many will be zero
      1. If you are looking at this post then we are measuring transfer functions from the blue flags to the green flags across the plant model.
      2. We will want to look at the TFs across the controller

 

  16461   Tue Nov 9 16:55:52 2021 Ian MacMillanSummaryComputer Scripts / ProgramsSUS Plant Plan for New Optics

[Ian, Tega]

We have moved c1sim computer from the test stand to the server rack in the office area. (see picture)

It is connected to the general campus network. Through the network switch at the top of the rack. This switch seeds the entire Martian network.

Test to show that I am not lying:

  1. you can ping it or ssh into it at
    controls@131.215.114.116
    Using the same password as before. Notice this is not going through the nodus network.
  2. It also has a different beginning of the IP addresses. Martian network IP addresses start with 191.168.113

c1sim is now as connected to the 40m network as my mom's 10-year-old laptop.

unfortunately, I have not been able to get the x2go client to connect to it. I will have to investigate further. It is nice to have access to the GUI of c1sim occasionally.

Attachment 1: IMG_8107.JPG
IMG_8107.JPG
  16462   Tue Nov 9 18:05:03 2021 Ian MacMillanSummaryComputer Scripts / ProgramsSUS Plant Plan for New Optics

[Ian, Tega]

Now that the computer is in its new rack I have copied over the filter two files that I will use in the plant and the controller from pianosa:/opt/rtcds/caltech/c1/chans to the docker system in c1sim:/home/controls/docker-cymac/chans. That is to say, C1SUP.txt -> X1SUP.txt and C1SUS.txt -> X1SUS_CP.txt, where we have updated the names of the plant and controller inside the txt files to match our testing system, e.g. ITMX -> OPT_PLANT in plant model and ITMX -> OPT_CTRL in the controller and the remaining optics (BS, ITMY, PRM, SRM) are stripped out of C1SUS.txt in order to make X1SUS_CP.txt. 

Once the filter files were copied over need to add them to the filters that are in my models to do this I run the commands:

$  cd docker-cymac
$  eval $(./env_cymac)
$  ./login_cymac
 #  cd /opt/rtcds/tst/x1/medm/x1sus_cp
 #  medm -x X1SUS_OPT_PLANT_TM_RESP.adl

see this post for more detail

Unfortunately, the graphics forwarding from the docker is not working and is giving the errors:

arg: X1SUS_OPT_PLANT_TM_RESP.adl

locateResource 'X1SUS_OPT_PLANT_TM_RESP.adl'

isNetworkRequest X1SUS_OPT_PLANT_TM_RESP.adl

canAccess('X1SUS_OPT_PLANT_TM_RESP.adl', 4) = 0

can directly access 'X1SUS_OPT_PLANT_TM_RESP.adl'

isNetworkRequest X1SUS_OPT_PLANT_TM_RESP.adl

locateResource(X1SUS_OPT_PLANT_TM_RESP.adl...) returning 1

Error: Can't open display:

This means that the easiest way to add the filters to the model is through the GUI that can be opened through X2go client. It is probably easiest to get that working. graphics forwarding from inside the docker is most likely very hard. 

unfortunately again x2go client won't connect even with updated IP and routing. It gives me the error: unable to execute: startkde. Going into the files on c1sim:/usr/bin and trying to start startkde by myself also did not work, telling me that there was no such thing even though it was right in front of me.

  16466   Mon Nov 15 15:12:28 2021 Ian MacMillanSummaryComputer Scripts / ProgramsSUS Plant Plan for New Optics

[Ian, Tega]

We are working on three fronts for the suspension plant model:

  1. Filters
    1. We now have the state-space matrices as given at the end of this post. From these matrices, we can derive transfer functions that can be used as filter inputs. For a procedure see HERE. We accomplish this using Matlab's built-in ss(A,B,C,D); function. then we make it discrete using c2d(sys, 1/f); this gives us our discrete system running at the right frequency. We can get the transfer functions of either of these systems using tf(sys);
    2. from there we can copy the transfer functions into our photon filters.  Tega is working on this right now.
  2. State-Space
    1. We have our matrices as listed at the end of this post. With those compiled into a discrete system in MatLab we can use the code Chris made called rtss.m to convert this system into a .c file and a .h file.
    2. from there we have moved those files under the userapps folder in the docker system. then we added a c-code block to our .mdl model for the plant and pointed it at the custom c file we made. See section 7.2 of T080135-v10
    3. We have done all this and this should implement a custom state-space function into our .mdl file. the downside of this is that to change our SS model we have to edit the matrices we can't edit this from an medm screen. We have to recompile every time.
  3. Python Check
    1. This python check is run by Raj and will take in the state-space matrices which are given then will take transfer functions along all inputs and outputs and will compare them to what we have from the CDS model.

 

Here are the State-space matrices:

A=\begin{bmatrix} 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ -\omega_x^2(1+i/Q_{x}) & -\gamma_x & \omega_xb & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ \frac{\omega_{\theta}^2}{l+b} & 0 & -\omega_{\theta}^2(1+i/Q_{\theta}) & -\gamma_{\theta} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & -\omega_{\phi}^2(1+i/Q_{\phi}) & -\gamma_{\phi} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & -\omega_{y}^2(1+i/Q_{y}) & -\gamma_{y}\end{bmatrix} 

  B=\begin{bmatrix} 0 & 0 & 0 & 0 \\ \frac{1}{m} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & \frac{R_m}{I_{\theta}} & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{R_m}{I_{\phi}} & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{m} \end{bmatrix}      C=\begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \end{bmatrix}      D=\begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix}

A few notes: If you want the values for these parameters see the .yml file or the State-space model file. I also haven't been able to find what exactly this s is in the matrices.

UPDATE [11/16/21 4:26pm]: I updated the matrices to make them more general and eliminate the "s" that I couldn't identify. 

The input vector will take the form:

\begin{bmatrix} x \\ \dot{x} \\ \theta \\ \dot{\theta} \\ \phi \\ \dot{\phi} \\ y \\ \dot{y} \end{bmatrix}

where x is the position, theta is the pitch, phi is the yaw, and y is the y-direction displacement

 

 

  16469   Tue Nov 16 17:29:49 2021 Ian MacMillanSummaryComputer Scripts / ProgramsSUS Plant Plan for New Optics

[Ian, Tega]

Updated A, B, C, D matrices for the state-space model to remove bugs in the previous estimate of the system dynamics. Updated the last post to represent the current matrixes.

We used MatLab to get the correct time-series filter coefficients in ZPK format and added them to the filters running in the TM_RESP filter matrix.

Get the pos-pos transfer function from the CDS model. Strangely, this seems to take a lot longer than anticipated to generate the transfer function, even though we are mainly probing the low-frequency behavior of the system.  

For example, a test that should be taking approximately 6 minutes is taking well over an hour to complete. This swept sine (results below) was on the low settings to get a fast answer and it looks bad. This is a VERY basic system it shouldn't be taking this long to complete a Swept sine TF.

 

Noticed that we need to run eval $(./env_cymac) every time we open a new terminal otherwise CDS doesn't work as expected. Since this has been the source of quite a few errors already, we have decided to put it in the startup .bashrc script.  

loc=$(pwd)
cd ${HOME}/docker-cymac/
eval $(./env_cymac)
cd ${loc}
Attachment 1: x_x_TF1.pdf
x_x_TF1.pdf
  16477   Thu Nov 18 20:00:43 2021 Ian MacMillanSummaryComputer Scripts / ProgramsSUS Plant Plan for New Optics

[Ian, Raj, Tega]

Here is the comparison between the results of Raj's python model and the transfer function measurement done on the plant model by Tega and me.

As You can see in the graphs there are a few small spots of disagreement but it doesn't look too serious. Next we will measure the signals flowing through the entire plant and controller.

For a nicer (and printable) version of these plots look in the zipped folder under Plots/Plant_TF_Individuals.pdf

Attachment 1: Final_Plant_Testing.zip
  16481   Wed Nov 24 11:02:23 2021 Ian MacMillanSummaryComputersQuantization Noise Calculation Summary

I added mpmath to the quantization noise code. mpmath allows me to specify the precision that I am using in calculations. I added this to both the IIR filters and the State-space models although I am only looking at the IIR filters here. I hope to look at the state-space model soon. 

Notebook Summary:

I also added a new notebook which you can find HERE. This notebook creates a signal by summing two sine waves and windowing them.

Then that signal is passed through our filter that has been limited to a specific precision. In our case, we pass the same signal through a number of filters at different precisions.

Next, we take the output from the filter with the highest precision, because this one should have the lowest quantization noise by a significant margin, and we subtract the outputs of the lower precision filters from it. In summary, we are subtracting a clean signal from a noisy signal; because the underlying signal is the same, when we subtract them the only thing that should be left is noise. and since this system is purely digital and theoretical the limiting noise should be quantization noise.

Now we have a time series of the noise for each precision level (except for our highest precision level but that is because we are defining it as noiseless). From here we take a power spectrum of the result and plot it.

After this, we can calculate a frequency-dependent SNR and plot it. I also calculated values for the SNR at the frequencies of our two inputs. 

This is the procedure taken in the notebook and the results are shown below.

Analysis of Results:

The first thing we can see is that the precision levels 256 and 128 bits are not shown on our graph. the 256-bit signal was our clean signal so it was defined to have no noise so it cant be plotted. The 128-bit signal should have some quantization noise but I checked the output array and it contained all zeros. after further investigation, I found that the quantization noise was so small that when the result was being handed over from mpmath to the general python code it was rounding those numbers to zero. To overcome this issue I would have to keep the array as a mpmath object the entire time. I don't think this is useful because matplotlib probably couldn't handle it and it would be easier to just rewrite the code in C. 

The next thing to notice is sort of a sanity check thing. In general, low precision filters yield higher noise than high precision. This is a good quick sanity check. However, this does not hold true at the low end. we can see that 16-bit actually has the highest noise for most of the range. Chris pointed out that at low precisions that quantization noise can become so large that it is no longer a linearly coupled noise source. He also noted that this is prone to happen for low precision coefficients with features far below the Nyquist frequency like I have here. This is one explanation that seems to explain the data especially because this ambiguity is happening at 16-bit and lower as he points out. 

Another thing that I must mention, even if it is just a note to future readers, is that quantization noise is input dependent. by changing the input signal I see different degrees of quantization noise.

Analysis of SNR:

One of the things we hoped to accomplish in the original plan was to play around with the input and see how the results changed. I mainly looked at how the amplitude of the input signal scaled the SNR of the output. Below I include a table of the results. These results were taken from the SNR calculated at the first peak (see the last code block in the notebook) with the amplitude of the given sine wave given at the top of each column. this amplitude was given to both of the two sine waves even though only the first one was reported. To see an example, currently, the notebook is set up for measurement of input amplitude 10.

  0.1 Amplitude of input 1 Amplitude 100 Amplitude 1000 Amplitude
4-bit SNR 5.06e5 5.07e5 5.07e5 5.07e5
8-bit SNR 5.08e5 5.08e5 5.08e5 5.08e5
16-bit SNR 2.57e6 8.39e6 3.94e6 1.27e6
32-bit SNR 7.20e17 6.31e17 1.311e18 1.86e18
64-bit SNR 6.0e32 1.28e32 1.06e32 2.42e32
128-bit SNR unknown unknown unknown unknown

As we can see from the table above the SNR does not seem to relate to the amplitude of the input. in multiple instances, the SNR dips or peaks in the middle of our amplitude range.

 

Attachment 1: PSD_IIR_all.pdf
PSD_IIR_all.pdf
  11389   Wed Jul 1 16:16:46 2015 IgnacioUpdateGeneralAccelerometers reinstalled for future huddle test

Today, I installed the Wilcoxon accelerometers in the table near the end of the mode cleaner. I only set three of them up instead of all six. They were set up just as Rana suggeted we should have them properly set up, i.e. cables being tighten up, and a box on top to prevent any airflow introducing any disturbances. We are planning on running the huddle test on these guys once the upgrade? to the interferometer is done.

The cables were tightly clamped to the table as shown below, I used a thick piece of shock absorbing rubber to do this.

 A small piece of thin rubber was used to hold each of the accelerometers tightly to the table in order not to damage them.

We had to borrow Megan's and Kate's piece of black foam in order to seal one of the sides properly, as the cable had to come out through somewhere. We didn't want to mess with drilling any holes into the box! 

There was a small crack even after using the foam. I sealed it up with duck tape.

The box isn't perfect, so there were multiple cracks along the bottom and top of it that could potentially allow for air to flow to the inside. Eric suggested that we should be super careful this time and do it right, so every crack was sealed up with ducktape.

 

Finally, we needed something heavy to be placed on top of the box to hold everything well. We used Rana's baby to accomplish this goal.

Just kidding! Rana's baby is too delicate for our purposes. A layman box of heavy cables was used instead. 

 

  11391   Sun Jul 5 18:14:13 2015 IgnacioUpdatePEMWilcoxon Accelerometer Huddle Test

Updated: On Thursday/Friday (sorry for late elog) I was messing with Eric's Wilcoxon 731A accelerometer huddle test data that was taken without the box and cables being adjusted properly. Anyways, I performed the three cornered hat analysis as he had done but I also performed a six cornered hat method as well instead of permuting around in pairs of three accelerometers. The following plots of the ASD's show the results,

It is interesting to note the improvement at low frequencies when six accelerometers are used instead of six while at higher frequencies we can clearly see how the results are worst than the three hat results.

I decided to take a mean of all six accelerometers measured ground signal as well as that for their computed selfnoises, this is plotted below,

 

Notice the obvious improvement along the entire frequency band of the measurements when all accelerometers are used in the data analysis.

I also performed some Wiener filtering of this data. There was an obvious improvement in the results,

The mean of the signals is also plotted below, just as I did with the cornered hat methods,

 

I also compared the mean self noise of the accelerometers against the manufacturers calculated selfnoise that Rana put up on Github. Both methods are compared against what the manufacturer claims,

As expected the measured noise curves of the Wilcoxon is not as good as what the manufactures stated. This should improve once we redo the huddle test with a better experimental setup. We have some pretty interesting results with the six cornered hat method at around 5 Hz, it is surprisingly accurate given how rough the calculations seemed to be.

I have attached my code for reference: code_accel.zipselfnoise_allsix.png

SEE attachments for better plots of the six accelerometers...

Attachment 5: code_accel.zip
Attachment 6: selfnoise_allsix_miso.png
selfnoise_allsix_miso.png
Attachment 8: selfnoise_allsix.png
selfnoise_allsix.png
  11407   Tue Jul 14 10:23:27 2015 IgnacioUpdateGeneralOptimal detector array placement thoughts

Over the past few days, I've been thinking about how to workout the details conerning Rana's request about a 'map' of the vicinity of the 40m interferometer. This map will take the positions of N randomly placed seismic sensors as well as the signals measured by each one of them and the calculated cross correlations between the sensors and between the sensors and the test mass of interest to give out a displacement vector with new sensor positions that are close to optimum for better seismic (and Newtonian) noise cancellation.

Now, I believe that much of the mathematical details have been already work out by Jenne in her thesis. She explains that the quantity of interest that we wish to minimize in order to find an optimal array is the following,

R = \sqrt{1-\frac{\vec{C}_{SN}^T C_{SS}^{-1}\vec{C}_{SN} }{C_{NN}}}

where  \vec{C}_{SN} is the cross-correlation vector between the seismic detectors and the seismic (or Newtonian) noise, C_{SS} is the cross-correlation matrix between the sensors and C_{NN} is the seismic (or Newtonian) noise variance. 

I looked at the paper that Jenne cited from which she obtained the above quantity and noted that it is a bit different as it contains an extra term inside the square root, it is given by

R' = \sqrt{1-\frac{\vec{C}_{SN}^T (C_{SS}^{-1}+C_{\Sigma\Sigma})\vec{C}_{SN} }{C_{NN}}}

where the new term, C_{\Sigma\Sigma} is the matrix describing the self noise of the sensors. I think Jenne set this term to zero since we can always perform a huddle test on our detectors and know the self noise, thus effectively subtracting it from the signals of interest that we use to calculate the other cross correlation quantities.

Anyways, the quantity R above is a function of the positions of the sensors. In order to apply it to our situation, I'm planning on:

     1) Performing the huddle tests on our sensors, redoing it for the accelerometers and then the seismometers (once the data aquisition system is working... )  

     2) Randomly (well not randomly, there are some assumptions we can make as to what might work best in terms of sensor placement) place the sensors around the interferometer. I'm planning on using all six Wilcoxon 731A accelerometers, the two Guralps and the STS seismometer (any more?).

     3) Measure the ground signals and use wiener filtering in order to cancel out their self noises.

     4) From the measured signals and their present positions we should be able to figure out where to move the sensors in order to optimize subtraction.

i have also been messing around with Jenne's code on seismic field simulations with the hopes of simulating a version of the seismic field around the 40m in order to understand the NN of the site a little better... maybe. While the data aquisition gets back to a working state, I'm planning on using my simulated NN curve as a way to play around with sensor optimization before its done experimentally.

i have as well been thinking and learning a little bit about source characterization through machine learning methods, specially using neural networks as Masha did back in her SURF project on 2012. I have also been looking at Support vector machines. The reasons why I have been looking at machine learning algorithms is because of the nature of the everchanging seismic field around the interferometer. Suppose we find a pretty good sensor array that we like. How do we make sure that this array is any good at some time t after it has been found? If the array mostly deals with the usual seismic background (quiet) of the site of interest, we could incorporate machine learning techniques in order to mitigate any of the more random disturbances that happen around the sites, like delivery trucks, earthquakes, etc.

  11423   Fri Jul 17 02:46:07 2015 IgnacioUpdateGeneralNew huddle test data for Wilcoxon 731A results

On Thursday, new huddle test data for the Wilcoxon 731A was aquired by Eric. 

The difference between this new data and the previous data, is:

1) We used three accelerometers instead of six this time around.

2) We used a foam box, and clamped cables on the experimental set up as shown in the previous elog, http://nodus.ligo.caltech.edu:8080/40m/11389

I have analyzed the new data. Here I present my results.

The following plot shows the ASD's for the three accelerometers raw outputs as well as their error signals computed using the three cornered hat method,

As before, I computed the mean for the output signals of the accelerometers above as well as their mean self noise to get the following plot

,

Now, below I compare the new results with the results that I got from the old data, 

Did the enclosure and cable clamping do much? Not really, according to the computed three hat results. Also, notice how much better, even if its a small improvement, we get from using six accelerometers and calculating their self noise by the six cornered hat method.

 

Now, I moved on to analyzing the same data with Wiener Filtering.

Here are again, the raw outputs, and the self noises of each individual accelerometer calculated using Wiener filtering,

The accelerometer in the Y direction is show a kind of funky signal at low frequncies. Why? Anyways, I calculated the mean of the above signals as I did for the three cornered hat method above to get the following, I also show the means of the signals computed with the old data using wiener filtering,

Is the enclosure really doing much? The Wiener filter that I applied to the huddle test old data gave me a much better, by an order of magnitude better self noise curve. Keep in mind that this was using SIX accelerometers, not THREE as we did this time. I want to REDO the huddle test for the WIlcoxon accelerometers using SIX accelerometers with the improved experimental setup to see what I get.

Finally, I compare the computed self noises above with what the manufacturer gives,

,

As I expected, the self noise using six accelerometers and Wiener filtering is the best I could work out. The three cornered hat method works out pretty well from 1 to 10 Hz, but the noise is just too much anywhere higher than 10 Hz. The enclosed, clamped, 3 accelerometer wiener filter result is an order of magnitude worse than the six accelerometer wiener filtered result, and two orders of magnitude worse than the three cornered hat method in the 1 to 10 Hz frequency band. 

As I stated, I think we must performed the huddle test with SIX accelerometers and see what kind of results we get.

Attachment 1: selfnoise_allthree_threehat_enclosed.png
selfnoise_allthree_threehat_enclosed.png
Attachment 2: selfnoise_3hat_enclosed_averages.png
selfnoise_3hat_enclosed_averages.png
Attachment 3: selfnoise_3hat_6hat_enc.png
selfnoise_3hat_6hat_enc.png
Attachment 4: miso_wiener_enclosedall.png
miso_wiener_enclosedall.png
Attachment 5: selfnoise_wiener_enclosed.png
selfnoise_wiener_enclosed.png
Attachment 6: compare_encl.png
compare_encl.png
  11424   Fri Jul 17 04:56:37 2015 IgnacioUpdateGeneralMCL Wiener filtering + FIR to IIR conversion using vectfit

 

We took data for the mode cleaner a while ago, June 30th I believe. This data contained signals from the six accelerometers and the three seismometers. In here I have only focused on the seimometer signals as witnesses in order to construct Wiener filters for each of the three seismometer signals (x,y,z) and for the combined seismometers signal. The following plot showing the ASD's shows the results,

 Wiener filtering works beautifully for the seismometers. Note that subtraction is best when we use all three seismometers as the witnesses in the Wiener filter calculation, as can be clearly seen in the first plot above.

Now, I used vectfit to conver the Wiener FIR filters for each seismometer to their IIR versions. The following are the bode plots for the IIR filters,

For the x-direction seismometer,

 

For the y-direction seismometer

,

 

And for the z-direction seismometer,

 The IIR filters were computed using 5 zeros and 5 poles using vectfit. That was the maximum number of poles that I could use wihtout running into trouble with matrices being almost singular in Matlab. I still need to figure out how to deal with this issue in more detail as fitting the y-seismometer was a bit problematic. I think having a greater number of poles will make the fitting a bit easier.

Attachment 1: Wiener_MCL_seismometers.png
Wiener_MCL_seismometers.png
Attachment 2: seisx_mag.png
seisx_mag.png
Attachment 3: seisx_mag.png
seisx_mag.png
Attachment 4: seisx_mag.png
seisx_mag.png
Attachment 5: seisx_phase.png
seisx_phase.png
Attachment 6: seisy_mag.png
seisy_mag.png
Attachment 7: seisy_phase.png
seisy_phase.png
Attachment 8: seisz_mag.png
seisz_mag.png
Attachment 9: seisz_phase.png
seisz_phase.png
  11425   Sat Jul 18 06:12:07 2015 IgnacioUpdateGeneralMCL Wiener filtering + FIR to IIR conversion using vectfit (Update)

(updateAfter Eric gave me feedback on my previous elog post, I went back and fixed some of the silly stuff I stated.

First of all, I have come to realized that it makes zero sense to plot the ASD's of the mode cleaner against the seismometer noise. These measurements are not only quite different, but elementary, they posess different units. I have focused my attention to the MCL being Wiener filtered with the three siesmometer signals. 

One of the major improvements that I make in the following analysis is,

1) Prefiltering; a band pass filter from 1 to 5 Hz, in order to emphasize subtraction of the bump shown in the figure below.

2) I have used vectfit exclusively in the 1 to ~5 Hz range, in order to model the FIR filter properly, as in, the kind of subtraction that we care about. Limiting myself to the 1 - 5 Hz range has allowed me to play freeley with the number of poles, hence being able to fit the FIIR filter properly with an IIR rational transfer function properly,

The resulting ASD's are shown below, in blue we show the raw MCL output, in blac the Wiener filter (FIR) result, and finally in black, the resultant data being filtered with the calculated IIR Wiener filter.

 

Now, in the following plots I show the IIR Wiener filters for each of the three seismometers, 

X Seismometer,

For the Y seismometer,

and for the Z seismometer,

 

The matlab code for this work is attached: code.zip

Attachment 1: Wiener_MCL_seismometers_iir.png
Wiener_MCL_seismometers_iir.png
Attachment 2: seisx_mag.png
seisx_mag.png
Attachment 3: seisx_mag.png
seisx_mag.png
Attachment 4: seisx_mag.png
seisx_mag.png
Attachment 5: seisx_ph.png
seisx_ph.png
Attachment 6: seisy_mag.png
seisy_mag.png
Attachment 7: seisy_mag.png
seisy_mag.png
Attachment 8: seisy_mag.png
seisy_mag.png
Attachment 9: seisy_ph.png
seisy_ph.png
Attachment 10: seisz_ph.png
seisz_ph.png
Attachment 11: seisz_ph.png
seisz_ph.png
Attachment 12: code.zip
Attachment 13: seisz_mag.png
seisz_mag.png
Attachment 14: seisz_mag.png
seisz_mag.png
Attachment 15: seisz_ph.png
seisz_ph.png
  11432   Tue Jul 21 05:17:09 2015 IgnacioUpdateGeneralMore clear accelerometer huddle tests results

I generated the following plots from the two sets of huddle test data we have for the accelerometers. 

Old data: 6 accelerometers, no cables clamped, no box, 55 mins worth of data.

New data: 3 accelerometers, cables clamped, foam box put on placed and completely sealed, 20 mins worth of data.

I made sure to use the same Impuse response time (6 sec) and sampling frequency (256 Hz), as well as every other parameter for the calculations.

  

 

Top left: The resultant self noise curve using the new data, there is definitely and improvement in the 0.5-2 Hz band. 

Top right: Resultant self noise using the old data, for the first set of three accelerometers.

Bottom left: Old data result for the remaining three accelerometers.

Bottom right: Old data result, using all six accelerometers as witnesses instead.

Attachment 1: new_data.png
new_data.png
Attachment 2: new_data.png
new_data.png
Attachment 3: old_data_1.png
old_data_1.png
Attachment 4: old_data_1.png
old_data_1.png
Attachment 5: old_data_2.png
old_data_2.png
  11442   Thu Jul 23 21:12:14 2015 IgnacioUpdateGeneralNewer accelerometer huddle test data + detrending

In the last post concerning the self noise of the accelerometers, I mentioned the differences between the two data sets I was playing with. In order to give a more concrete analysis of the accelerometers self noise, we came to the conclusion that data taking time should be the same.

The plots below show the analysis for the following two datasets:

Old Data:  6 accelerometers, no cables clamped, no box, 55 mins worth of data.

Newer data: 3 accelerometers, cables clamped, foam box put on placed and completely sealed, 57.66 mins worth of data, (we had 20 mins of data in the previous data set).

Filter parameters were kept the same in all calculations, the only change that was added to the analysis was the detrending of the signals using the detrend function on Matlab, this improved the results as the plots show. I also plotted the error bars for the Wiener filtered result for reference as well as the manufactures claimed self noise.

   

After detrending the data and taking a longer dataset we can see the improvement brought upon by the foam box in the low frequency band of 0.5 - 10 Hz, as shown in the bottom left plot. There is a lot of noise that needs to be cancelled out from 10 Hz and on, which brings to our attention the plot on the bottom right corner. This plot uses the old data but uses all six accelerometers as witnesses, it also improved overall after having detrended the data, but what is peculiar about this plot is the fact that it manages to mitigate the higher frequency 10 - ~100 Hz band noise. 

 

Attachment 1: old_data_1.png
old_data_1.png
Attachment 2: old_data_2.png
old_data_2.png
Attachment 3: new_data.png
new_data.png
Attachment 4: six_old_data.png
six_old_data.png
  11446   Fri Jul 24 23:08:53 2015 IgnacioUpdatePEMMC1 accelerators moved for future huddle test

I have moved the MC1 accelerators and cable to the huddle test set up, in order to see how a six witness huddle test with the improved set up will do. 

Here is a picture of the accelerometer set up,

Our motivation for doing this is to see if more witness signals used in the Wiener filter really does indeed improve subtraction, as it was seen from previous huddle results, specially in the region above 10 Hz.

  11449   Tue Jul 28 05:00:03 2015 IgnacioUpdateGeneralSeismometer cans

I've have been talking a little bit with Steve about the seismometer enclosures.

We want to improve on the current stainless steel cans that cover the two Guralps at the end of the arms. In order to do this, we want to cover the interior of the cans with copper foil to improve the thermal conductivity of the enclosure to better control the temperature inside it. Ideally, we would want to copper plate the cans, but cost and difficulty has been an issue.

I have done some rough calculations and it seems that we need a copper layer of thickness being about a third that of the stainless steel can. This happens to be around 0.5-0.6 mm since we have 16 gauge (~1.6 mm) stainless steel cans. 

After wrapping the cans interior with copper, we will insulate them with foam in order to improve its thermal inertia. We want to probably use the same foam that Megan has been using for her seismometer enclosure. I have yet to think about a heater, but something similar to Megans resistor thing would work only smaller. I would be placed inside the can, right on the center of its bottom in order to ditribute heat evenly.

 

  11450   Tue Jul 28 09:31:35 2015 IgnacioUpdateGeneralNewest accelerometer huddle test

I downloaded new accelerometer huddle test data from last night and analyzed it. This new data set is ~55 mins and uses the same Wiener filter parameters as previous huddle test analysis, the main difference being six accelerometers used in the Wiener filter and the improved experimental set up.

After computing the ASD for the self noise for each of the six accelerometers, (being witnessed by the remaining five), we get,

Now computing the mean of the above signals and the corresponding error bars gives the following result,

Comparing to prevoius huddle tests, I can note two trends on the Wiener subtraction:

1) When using six accelerometers, the subtraction above ~8 Hz drastically improves.

2) When using three accelerometers, there is better cancellation in the 1-5 Hz region, see http://nodus.ligo.caltech.edu:8080/40m/11442. This might have to do with how much more closer the accelerometers are to each other? 

Attachment 1: selfnoise_allsix_miso_newest.png
selfnoise_allsix_miso_newest.png
  11457   Wed Jul 29 10:34:42 2015 IgnacioSummaryLSCCoherence of arms and seismometers

Jessica and I took 45 mins  (GPS times from 1122099200 to 1122101950) worth of data from the following channels:

C1:IOO-MC_L_DQ (mode cleaner)
C1:LSC-XARM_IN1_DQ (X arm length)
C1:LSC-YARM_IN1_DQ (Y arm length)

and for the STS, GUR1, and GUR2 seismometer signals.

The PSD for MCL and the arm length signals is shown below,

I looked at the coherence between the arm length and each of the three seismometers, plot overload incoming below,

For the coherence between STS and XARM and YARM,

  

For GUR1,

  

Finally for GUR2,

 

A few remarks:

1) From the coherence plots, we can see that the arm length signals are coherent with the seismometer signals the most from 0.5 - 50 Hz. This is most evident in the coherence with STS. I think subtraction will be most useful in this range. This agrees with what we see in the PSD of the arm length signals, the magnitude of the PSD starts increasing from 1 Hz and reaches a maximum at about 30 Hz. This is indicative of which frequencies most of the noise is present.

2) Eric did not remember which of  GUR1 and GUR2 corresponded to the ends of XARM and YARM. So, I went to the end of XARM, and jumped for a couple seconds to disturb whatever Gurald was in there. Using dataviewer I determined it was GUR1. Anyways, my point is, why is GUR1 less coherent with both arms and not just XARM?  Since it is at the end of XARM, I was expecting GUR1 to be more coherent with XARM. Is it because, though different arms, the PSD's of both arms are roughly the same? 

3) Similarly, GUR2 shows about the same levels of coherence for both arms, but it is more coherent. Is GUR2 noisier because of its location?

Code: ARMS_COH.m.zip

Attachment 1: PSD_ARMS_MCL.png
PSD_ARMS_MCL.png
Attachment 2: XARM_STS_COH.png
XARM_STS_COH.png
Attachment 3: YARM_STS_COH.png
YARM_STS_COH.png
Attachment 4: XARM_GUR1_COH.png
XARM_GUR1_COH.png
Attachment 5: YARM_GUR1_COH.png
YARM_GUR1_COH.png
Attachment 6: XARM_GUR2_COH.png
XARM_GUR2_COH.png
Attachment 7: YARM_GUR2_COH.png
YARM_GUR2_COH.png
Attachment 8: ARMS_COH.m.zip
  11460   Wed Jul 29 17:51:56 2015 IgnacioUpdatePEMAccelerators moved back to MC1 and MC2

We are done taking accelerator huddle test data. So I moved back all six accelerometers and cables to MC1 and MC2. I also relabel each of the accelerometers properly since the labels on them were confusing.

QED

 

Attachment 1: FullSizeRender.jpg
FullSizeRender.jpg
Attachment 2: FullSizeRender_2.jpg
FullSizeRender_2.jpg
  11462   Thu Jul 30 02:06:20 2015 IgnacioUpdateIOOMC2 <-> MCL Actuator TF fitted

Eric downloaded MC2 to MCL transfer function data (H) as well as its inverse, MCL to MC2 (Hinv). He also downloaded new MCL and MC2 data.

I used vectfit to fit the MC2 to MCL transfer function, 

The ZPK parameters for this fit were,

Zeros              1278.36719876674 + 0.00000000000000i
                   -100.753249679343 + 0.00000000000000i
                   -18.6014192997845 + 13.0294910760217i
                   -18.6014192997845 - 13.0294910760217i

Poles              -1.11035771175328 + 7.03549674098987i
                   -1.11035771175328 - 7.03549674098987i
                   -18.8655320274072 + 0.00000000000000i
                   -690.294337433234 + 0.00000000000000i

Gain               0.00207206036014220

Using the above vectfit model, I filtered the raw MC2 signal to get 'MCL'. The PSD's of the raw MCL data and the filtered MC2 result is shown below,

The lack of accuracy of the transfer function at replicating MCL at frequencies lower than 0.7Hz is expected, the vectfit model I generated fails to follow accurately the raw transfer function data. My question: Does it matter? My guess: Probably not. In order to mitigate seismic noise from the mode cleaner we are mainly concerened with the 1-3 Hz region.

I also used vectfit to fit the transfer function for MCL to MC2,

This one was harder to fit accurately for some reason, I could do it with four pairs of zeros and poles but it took some preweighting.

The ZPK parameters for the above fit were, 

Zeros              0.173068278283995 + 0.00000000000000i
                   0.995140531040529 + 0.0268079821980457i
                   0.995140531040529 - 0.0268079821980457i
                   0.894476816129099 + 0.00000000000000i

Poles              -19.9566906920707 + 18.0649464375308i
                   -19.9566906920707 - 18.0649464375308i
                   -109.275971483008 + 0.00000000000000i
                   -1791.88947801703 + 0.00000000000000i

Gain               1237.46417532120

Similarly, using this ZPK model, I filtered the MCL signal to get 'MC2'. I plotted the PSD for the MC2 signal and the filtered MCL to get,

Again, the lack of accuracy of the filtered MC2 at replicating MCL below 0.7 Hz and above 12 Hz is due to the inverse transfer function failing to converge in these ranges.

Attachment 1: TF_BODE.png
TF_BODE.png
Attachment 2: MC2_2_MCL.png
MC2_2_MCL.png
Attachment 3: TF_INV_BODE.png
TF_INV_BODE.png
Attachment 4: MCL_2_MC2.png
MCL_2_MC2.png
  11472   Thu Jul 30 19:12:52 2015 IgnacioUpdateIOOYAW and PIT WFS Wiener filtering

Rana pointed out that another way to mitigate seismic motion at in the mode cleaner would be to look at the YAW and PITCH output  channels of the WFS sensors that control the angular alignment of the mode cleaner. 

I downloaded 45 mins of data from the following two channels:

C1:IOO-WFS1_YAW_OUT_DQ

C1:IOO-WFS1_PIT_OUT_DQ

And did some quick offline Wiener filtering with no preweighting, the results are shown in the PSD's below,

and

I'm quite surprised at the Wiener subtraction obtained for the YAW signal, it required no preweighting and there is about an order of magnitude improvement in our region of interest, 1-3 Hz. The PIT channel didn't do so bad either.

 

Attachment 1: YAW.png
YAW.png
Attachment 2: PIT.png
PIT.png
  11485   Thu Aug 6 21:03:45 2015 IgnacioHowToWienerFilteringHow to do online static IIR Wiener filtering

In order to do online static IIR Wiener filtering one needs to do the following,

1) Get data for FIR Wiener filter witnesses and target.

2) Measure the transfer function (needs to be highly coherent, ~ 0.95 for all points) from the actuactor to the control signal of interest (ie. from MC2_OUT to MC_L).

3) Invert the actuator transfer function.

4) Use Vectfit or (LISO) to find a ZPK model for the actuator transfer and inverse transfer functions.

5) Prefilter your witness data with the actuator transfer function, to take into account the actuator to control transfer function when designing the offline Wiener FIR filter.

6) Calculate the frequency response for each witness from the FIR coefficients.

7) Vectfit the frequency reponses to a ZPK model, this is the FIR to IIR Wiener conversion step.

8) Now, either, divide the IIR transfer function by the actuator transfer function or more preferably, multiply by the inverse transfer function.

9) Use Quack to make SOS model of the IIR Wiener filter and inverse transfer function product that goes into foton for online implementation.

10) Load it into the system.

The block diagram below summarizes the steps above.

Attachment 1: iir.png
iir.png
  11488   Mon Aug 10 22:18:19 2015 IgnacioUpdateIOOReady to do some online mode cleaner subtraction

I'm attaching a SISO IIR Wiener filter here for reference purposes that will go online either tonight or tomorrow evening. This is a first test to convince myself that I can get this to work, MISO IIR filters are close to being ready and will soon be employed. 

This Wiener filter uses the STS-X channel as a witness and MCL as target. The bode plot for the filter is shown below,

The performance of the FIR and IIR Wiener filters and the ammount of subtraction achive for MCL is shown below,

 

Output from quack to be loaded with foton: filter.zip

K bye.

Attachment 1: stsx.png
stsx.png
Attachment 2: performance.png
performance.png
Attachment 3: filter.zip
  11492   Tue Aug 11 11:30:19 2015 IgnacioUpdateIOOSISO (T240-X) FF of MCL

Last night we finally got some online subtraction going. The filter used is described in the post this eLOG is @eLOG 11488

The results were as follow:

The filter worked as expected when subtracting noise out of MCL,

There is about a factor of 6 subtraction at the ~3Hz resonant peak. The static IIR filter predicted a factor of 6-7 subtraction of this peak as well.

The 1.2 Hz resenonant feature improved by a factor of 3. This should improve quite drastically when I implement the y-channel of the T240 seismo.

There is some high frequency noise being injected, not very noticeable, but present. 

We then took a look at the power in the MC when the filter was on,

The power being transmitted in the cavity was not as stable as with the feedforward on. We believe that the filter is not at fault for this as Eric mentioned to me that the MC2 actuator lacked some sort of compensation that I need to understand a bit better.

YARM was then locked when the filter was on and we took a look at how it was doing. There was stationary sound arising from the locking of the YARM, leading us to believe that the filter might have injected some noise in the signal. IT DID.

The filter injected nasty high frequency noise at YARM from 11 Hz and on. This is to be expected since the filter did not roll off to zero at high frequencies. Implementing a 1/f rolloff should mitigate some of the injected noise.

 Also, as one can see above, subtraction by around a factor of 2 or so, was induced by the mode cleaner feedforward subtraction.

Attachment 1: MCL.png
MCL.png
Attachment 2: MCTRANS.png
MCTRANS.png
Attachment 3: YARM.png
YARM.png
ELOG V3.1.3-