IIR Filtering Ideas

This post is going to be in pseudo-formal speak since I'm fleshing out ideas for my dissertation. Although that's tough without adding any equations.

The problem of constructing the optimal IIR controller for a given closed-loop plant and disturbance model can be solved as an LQR problem with a particular state-space system. However, the condition that the plant transfer matrix is commutable with the filter turns out to be overly restrictive. Suppose the transfer matrix can be factored into the product of minimum and non-minimnum phase matrices. We can define a new filter which is the product of the minimum phase component and the optimal filter F, this leaving the non-minimum phase component to be compensated. The requirement is now that this non-minimum phase component is commutable, which happens if its equal to a scalar transfer function times the identify matrix. If this is satisfied, we can perform the LQR problem to identify the controller, then multiply it by the inverse of the minimum phase component to recover the actual optimal filter that is implemented in the software.

For the adaptive optics experiment, the lack of significant DM dynamics simplifies the problem. If the non-minimum phase transfer matrix consists solely of n-step delays on the diagonal, then the optimal IIR filter is simply the n-step Kalman predictor for the disturbance model. The disturbance model is itself in innovations form, thus the Kalman predictor can be constructed directly from the state-space model generated by the subspace identification algorithm.

...some math showing how this is done....

The result is a filter which predicts the disturbance wavefronts n-steps ahead on the basis of the current wavefront measurement.

Beautiful.

Overly Ambitious

Now that the FIR filter is basically working for the single channel experiment, I think its a good time to set some targets for the next four weeks. The clear next step is to start working on a multichannel version using more modes. Theoretically there isn't that much difference here, but I suspect there will be practical problems with actuator saturation and other shittyness that will slow things down.

I think a reasonable goal is a 10 channel adaptive/optimal filter in 4 weeks. Here are some things that will need to happen:

1. Characterize closed-loop transfer matrix. How similar are the diagonal terms? How significant are the off diagonals? What's the best scalar transfer function approximation?

2. Get the adaptive controller working given an identified or ideal transfer matrix.

3. Write a script to compute the optimal multichannel FIR and IIR filters.

4. Get target camera working to compute Strehl ratios.

5. Find/write/steal a simulink block that can implement multichannel transfer functions. The pole at 1 in the pure integrator I'm using now might (will) lead to saturation when more modes are controlled.

6. Questions: How to the PSD's of the output channels compare? How does changing the number of modes alter steady-state performance?

A lot of this will require running many experiments or simulations, so there should be plenty of down time to pursue some theoretical stuff for the SISO case. Basically, I'd like to write m-files that do the following:

1. Compute the optimal IIR filter. How does performance compare to the FIR case?

2. Compute the optimal FIR gains using and RLS array algorithm.

3. Compute the optimal FIR gains using an RLS lattice filter.

4. Do something about implementing a Laguerre filter.

SInce these will all be m-files, I don't expect to implement these in the actual experiment right away. Mainly, I want to get some idea of the theoretical performance benefits by using increasingly complicated methods. Obviously I want to do all of this for the multichannel case one day, but the details of that are so complicated my head might explode first. We'll see.

Phucking Transpose

Yes, a freaking missing apostrophe was responsible for nearly a weeks delay. To try to narrow down the problem with the impulse response filter calculation, I was trying to match the results with the data driven m-file using made up disturbance models comprised of random matrices. It turned out that both filters were the same as long as the A matrix was diagonal. The only place in the code where this mattered were locations where I needed A transpose. The function call to Matlab's dlyap function was the freaking culprit. That piece of shit is one of those functions that's screwed me in the past, and of course the one time I wasn't careful it bits me in the ass.

Anyway, with the fix the impulse response method now returns basically the same FIR filter coefficients as using data. I also realized that the results from the adaptive loop I posted yesterday were crap. Somehow I was using the wrong model for the closed loop plant, chalk that up to shitty variable names. Here are the proper PSD's and modal outputs


My advisor thinks these results are stellar. Amazingly the PSD with the optimal FIR filter is pretty similar to the AO PSD, hopefully showing that my idea of dividing out the part of the disturbance cancelled by the classical loop is correct. In this case both the adaptive and fixed gain FIR filter are using 4 taps, so you might ask why the AO loop does better than the "optimal" filter at certain frequencies. The reason is that the adaptive loop can compensate somewhat for modeling error.

Right now I'm trying to run things without the shitty DM pause, which speeds things up to around 20Hz. I suspect the results won't be so peachy, but the time savings would be huge (8 minutes vs 40 for 10000 frames), and I wouldn't have to spend so much time quality time with youtube waiting for my experiments to run.

Also, I'd like to look at how the number of filter coefficients changes the steady-state performance, although I don't think adding many more taps will make much difference. Also I'd be cute to have the Strehl ratio performance to look at too.

Peace out.

Case of the Mondays

Here are the first results from using the "optimal" FIR filter, computed using the data-drive approach. Surprisingly it performs pretty well compared to the AO loop, the fact that its not spitting out absolute crap is a small miracle.



Here's a sampling of the modal output


I'm not sure what's going on with the AO loop. Looking at the commands with the adaptive loop closed shows lots of lower-end saturation going on about 1/3 of the time, so something is probably screwed up somewhere in the experiment. Hard to say at the moment since I'm running this remotely from home. The frustrating part is that it takes to friggin long to run an experiment, around 40 minutes for 10000 frames, that its easy to distract my already OCD mindset. I'm going to have to start doing this in simulations first.

This is good for a first step, but there are still some outstanding questions I'd like to look at this week. The first few deal with this optimal FIR filter calculation:

1. Determine what's really causing the difference between the data and impulse response drive methods to finding the optimal gains.

2. What's the real disturbance model that should be used in the calculation? What's the difference between computing it and identifying it from i/o data? Should it be SISO or MISO?

3. How does the filter order affect performance?

4. Write a script to determine the optimal IIR filter by solving an LQR problem.

Also, all of this stuff so far has been for the first focus mode. Sooner or later I'm going to have to do everything over again al MIMO, so I'd be nice to have some heads up if there are potential problems in the road. The first step is to look at the transfer matrix for multiple modes with the classical loop closed. Everything depends on this being diagonal with the same SISO tf on the diagonals. If this doesn't hold to a reasonable extent then there could be serious limitations. With that in mind:

1. How similar are the diagonal transfer functions for each mode? Models identified with significant saturation are garbage.

2. If they're all of the same form, but with a different gain, can the transfer matrix be factored into a single transfer function times a static gain matrix? If so, can this matrix just be incorporated into the poke matrix?

3. What's the difference between doing a MIMO subspace ID and multiple SISO id's?

4. Does the simulation even have enough accuracy to identify the model for multiple channels?

All this will be much faster if I just suck it up and do it in silico first.

8.16.10

Still having trouble getting decent results with the FIR filter. I still can't get the impulse response method to agree with what the least-squares solution spits out, even when I include the noise covariance in the state-space model. I might try a third method, there the filter coefficients are spit out from a finite-time LQR problem.

One thing I noticed is that the disturbance model I use isn't exactly the state-space system identified directly from the disturbances. Instead, its the part of the disturbance left over after going through the classical control loop. After untangling the block diagram, you end up dividing the original disturbance system by some transfer function involving the plant model. The problem I have is that this transfer function might not be exactly minimum phase, so you can end up with an unstable disturbance model to put into your FIR calculation.

I think its time to visit my advisor.

8.10.10 [2]

I feel like I'm loosing the script here, so a quick sitrep on what's going on: I currently have 2 different methods written to calculate the optimal FIR disturbance rejection filter. One uses state-space models of the disturbance and plant to generate impulse response sequences. These can be used to form a Weiner-Hopf problem and solved for the optimal coefficients a la earlier work we did on jitter control. The second method uses the models to simulate data directly. Given enough samples, the solution to another (similar) linear equation yields the coefficients.

Theoretically, both of these should give similar results (I think). But so far no luck. Here's a comparison of the output PSD with filters calculated using each method compared to using no filter (F=1).



I used the actual plant disturbance models I identified from the experiment. Using the data approach works pretty well, although its pretty cumbersome, and would be stupid with multiple modes. The method using the impulse responses, however, is just crap, clearly making things worse.

Confusingly, both methods crap out the identical filter with less complicated disturbance models. I think the problem is that the input noise covariance matrix isn't really accounted for in the state-space model of the disturbance. It comes into play in the data driven case since I have to use it to generate the input, but it doesn't show up directly in the impulse response as the moment. It should be easily, however, to incorporate it into the state-space model by multiplying the input matrix.

Thats the plan for the afternoon, as soon as I finish blogging here in a coffee shop.

8.10.10

After doing a system ID on the new disturbance model, I was finally able to get my "optimal" FIR filter working. Or rather, functional since in the experiment I get saturation and in the simulation I get crap. Right now I'm trying to write a similar script that finds the optimal coefficients by using simulated data directly in the least-squares problem, instead of just impulse response terms. We'll see if I get the same results as with the other method.

I can already see the next stop on this pain train. The adaptive filter basically solves this least-squares problem recursively using either LMS or RLS filter, so the obvious next step is to write my own adaptive code. The only piece missing is my ignorance about writing S-functions for Simulink, but there are hints that I might have to learn that eventually anyway.

One Thing Leads To Another

I think my plan to do a subspace ID on the actual disturbances as measured by the WFS is a good one, but it turned out that the code we have to do it requires the measurements to be on a square grid. Stupidly, I've been using a rectangular 19x22 measurement area on the WFS, with the DM area (found by adding all the influence functions) shifted to the right by around 5 columns. Putting disturbances from a state space model on this odd configuration required scaling the image, shifting it in an attempt to align it with the DM area, an only then performing a least-squares fit with the poke matrix. Naturally the results were shit, as shown in the video in 6.18.10.

My slope calculation code is ludicrously cumbersome, so it took a day to properly rewrite it to use a 19x19 grid (actually a down-sampled 38x38 grid) and make sure it was bug-free. Realigning the beam and removing tilt so that the area for both DMs was centered took another. My advisor commented that he doesn't know how I keep everything straight in my head with so much shit going on. I assured him I have no idea what I'm doing.

Anyway, its good to periodically tinker with the experiment and realign everything anyway to keep my monkey skills sharp. Today I put some state-space generated disturbances on the DM, which now only requires rescaling the state-space output to a slightly larger grid. The resulting phases look much better. Here's a comparison between (L to R) the state-space output (1 phase screen model on a 17x17 grid), the DM phase predicted using the poke matrix (on a 19x19 grid), and the actual phase measured by the WFS (with the bias removed of course).



If you squint you can actually see the similarity in the phases as it flows across the aperture. Either way it much cleaner than the random flatulence I was getting before.

Tomorrow I'll look at how this affects the performance of the adaptive loop. Now that the grid is square I can also go ahead with my original plan and compute an optimal FIR filter.

7.31.10

Ahh crap its August.

Finished some code to calculate the optimal FIR filter (in the SISO case) in the disturbance rejection problem by solving the Weiner-Hopf equations. I have no clue if it works, but it seemed to give non-bullshit answers at the few, simple problems I threw at it. To actually use it in the experiment I have to harness the disturbance state-space model, but because of all the reshaping and massaging required to get the disturbance commands I'm not sure the original model would be valid. Right now I'm working on identifying a new model based on the measured wavefronts, so the basic procedure retardedly involves pulling a subspace ID twice. At that point I'm probably just getting complete nonsense, but its worth a try.

7.26.10

I was able to get the experiment running around 20Hz with Simulink by removing the pause between DM commands. This is great, except it now looks like DM dynamics are definitely screwing things up.

Here's a plot comparing the PSD's from 10000 frames using the same disturbance input and several different plant models.


Clearly using the original plant model in the AO loop, which contains a single delay and PI controller, produces crap; its worse than using just the classical controller. From my ARX experiments with the DM, I found that running at full speed seems to add an additional delay, and sure enough multiplying this ideal plant model with a delay works much better. Unsurprisingly the best results come with identifying a plant model first.

None of these results are as good as what I get with the artificial pause between DM commands, but the time difference is significant (20Hz vs. 4 Hz). I get the feeling this is a fact I'll just have to live with for the time being.

Can I graduate now?

Finally finally have the adaptive loop working in the experiment, at least with a single mode. Turns out that the effect on the overall wavefront norm is difficult to see for this disturbance model unless you know its there. Its more clear if you just look the coefficient of the mode being controlled.

But the best way to observe the effectiveness is by comparing the PSDs of the modal coefficient. In the uncontrolled case (classical PI controller only), you can clearly see some color resulting from the disturbance input. In simulation the adaptive loop flattens this out somewhat:



To my surprise, the adaptive loop does an even better job whitening the PSD in the experiment:



Another thing to notice is that I compared the performance using the ideal and identified plants. The adaptive loop uses a model of the closed (classical) loop to estimate the disturbance input. In general we assume that the ideal plant is just an integral controller and a unit delay, since the phase reconstructor is chosen to be the pseudo-inverese of the modal poke matrix. To verify this I also identified a plant using n4sid and a few thousand samples of input/output data, and found the resulting transfer functions very identical. This is nice to know since the plant actually contains some dirty nonlinearities like saturation and rounding, so it looks like those aren't significant for now.

This is all for a single mode. The requirements on the plant in the multiply mode case are more stringent (i.e. a diagonal transfer matrix). Something is also causing this to run quite a bit slower than from an m-file, so that will take some coffee consumption to figure out as well.

7.12.10 [2]

Wow, 2 posts within 24 hours. This is what happens when you can't nail down consistent sleep patterns.

Advisor was intrigued by the comparison of the velocity estimates (see plot in previous post). Particularly, both lines have almost the same shape, and appear to only differ by the average slope, ie the velocity. I proposed that this is because the poke matrix doesn't account for the change in the beam size that's roughly proportional to the norm of the command vector. Relatively small perturbations are used to estimate the poke matrix, so the beam diameter is actually smaller for general random commands. A particular phase profile enters and exits the DM surface in the same time for either case, so if the actual beam diameter is smaller it means that the phase traverses fewer pixels for a given time period, ergo resulting in a lower velocity estimate.

I have no idea if that makes any sense since its around 2AM. In any case, its not clear there's anything I can really do about it other than estimate some correcting scale factor and apply that to every command sequence. Its also not clear if any of this velocity estimation stuff will make it in a paper or my dissertation, so I'm not sure its worth devoting my entire life to something that's basically a sideshow to the main event.

Anyway, tomorrow (today), I'd like to ignore this discrepancy for now and look at applying phases at twice to rate to see if I can measure twice the velocity. I plan on doing this by using the same sequence of phases from the SS model, but just applying every other command.

I'm also making some (theoretical) progress on implementing an optimal FIR filter in the actual AO experiment. I still have to think about what it means to do the calculation in the multi-channel case.

7.12.10

More stuff on the velocity estimation. I managed to run some disturbances on the experiments that originated from a state space model. Surprisingly, you can actually distinguish something that looks like "flow" in the resulting reconstructed phase measurements. As a bonus, the velocity seems to be relatively constant.

The velocity estimate from these measurements are less than what's predicted by putting the estimated phases (using the commands and the phase poke matrix) through the estimator, but just the fact that there's anything recognizable is a plus. I did, however, have to mask out only the center of the WFS image corresponding to the active region of the DM, everything outside of this is just distortion. Maybe something like this should be done in the AO loop as well.



Speaking of which, I really want to focus on getting the AO loop working in the experiment this week. So far I everything runs, but I haven't seen any improvement in the Strehl in either the Simulink experiment or simulation. A few things to try:

1. Compare predicted and actual disturbance measurements (w/ and w/o bias). This should pin down if the internal plant model is accurate. It should be after running n4sid on sample data. Theoretically, I think the MSE between these should converge at something like an exponential rate after the adaptive loop is closed.

2. Try different disturbance sources. Maybe the current SS system is just too close to white to be useful. Maybe try a simple FIR filter or different amplitudes.

3. Compute the optimal IIR and FIR filter using the known disturbance model and see if that makes any difference.

Getting this working, especially #3, is important. The stuff with the velocity and new SLM is just icing at the moment.

7.8.10

Things to do today and tomorrow:

- Rewrite correlation code to handle phase data on rectangular grids

- Map commands from a SS model to DM commands, and compare velocity estimates for SS model, output from SS model, projected DM surface profile w/ and w/o bias. Is there any correspondence at all between the DM surface velocity and the SS estimate?

- Moar subspace ID stuff. Calculating oblique projections using LQ factorization.

- Screw around with SLM more.

Fun.

7.7.10

Still progressing on this correlation/velocity estimate stuff. Clearly from the previous post, estimating the velocity at each stop is produces inconsistent results. I think the problem is that when the speed is relatively slow (<1 px/frame), consecutive frames look very similar, with only a few edge pixels changing. Thus there just isn't enough movement in the first few delays to show clear movement of the peak of the correlation image, maybe explaining why it takes a few delays for the estimated velocity to settle down to a reasonable number.

But based on how consistently the peak moves in that video I decided to just track its position as a function of the delay, instead of producing an estimate each time. Luckily, the peak position is very linear. Because the peak should move with the same velocity as the phase profile, I can do a linear regression and take the resulting slope as the velocity. The resulting estimate is close to the average of estimating a velocity each delay, but its more justified looking at a plot of the peak position vs delay.

This works equally well using the covariance matrix of a state space model to calculate the correlation image for each time lag. To make it even better I managed to vectorize the calculation so that Matlab doesn't yack all over the nested loops. It runs around 10x faster than before.

I'm not sure where all this is going. It seems to be working pretty well, but I'm not sure we'll be able to squeeze a paper out of it. Maybe if I manage to get results of flow using the WFS and the actual experiment that would be more interesting.

6.25.10

As I mentioned, one problem with the velocity estimation method I'be been playing around with is that the estimates aren't constant for autocorrelations over a different number of delays, even when only one phase screen is used to generate the data. The question is, is this a result of the method, or is something phucked up in the simulation data itself?

We needed a second set of data to test things out on, and I remembered that we got this CD with "challenge" data from a conference we attended a few months ago. Basically, the data are frames from OPD data of turbulence over some flat plate, so its not exactly the same as general AO turbulence, but close enough to validate the velocity estimation algorithm. The data's pretty dense: each frame has phases on a 41x41 grid, and there are around 15000 frames in each file. Unbelievably though, its all stored uncompressed in a fucking 1GB text file.

It took me 2 days to figure out how to load just enough data in from this shit pipe to be useful. Right now I'm taking every third frame or so, and only using every other grid point. Making Matlab do this without reading every line in the text file took more than a few hits of caffeine.

Just looking at the autocorrelations of the frames (instead of calculating a SS model first), I still see the problem if varying velocity estimates as a function of the delay. Calculating the velocity from the data itself is something other people have done successfully, so either they're all full of shit or I'm just missing some detail. I suspect there's some pre-processing of the data I have to do to get more consistent estimates.



Here's a cool video of the spatial autocorrelation over 7500 frames for a varying number of delays. The peak should move with the same velocity as as they layer, at least until enough delays are used so the frames are essentially uncorrelated.



From the vid it looks like its moving with a constant velocity, but its hard to know exactly.

I need a drink.

6.18.10

This is "inter-session" week, the time of bliss on campus between finals and the start of the summer quarter when undergrads have gone on internships and professors have gone on vacation, and only pale grad students and vagrants are out and about. Its basically the only time its possible to find an empty seat in the campus coffee shop, which staffed by dour students clearly disappointed they couldn't find something better to do for the summer.

You'd think with finals over that this would be a great time to for me to actually get some work done, but no, I had to move to another apartment just as it was starting. After a few days unpacking and waiting for utility people I'm finally starting to get back into the swing of things.

The velocity/correlation analysis I've been working on still shows some promise. I've tried it on "real" simulated data and it is still able to produce discernible peaks from the state-space model of the disturbance. The problem is that the resulting velocity estimate varies depending on the number of delays in the correlation, essentially implying that the velocity isn't constant across the frame. These simulation are supposedly using a single phase screen, but since they're essentially black boxes for us who really knows what's happening. Its hard to know then if my method is return accurate velocity estimates.

Right now the idea is simmering on the back burner until I can get more state-space models. Another option is to use this "challenge data" from a conference we went to that is apparently the turbulence from flow over some kind of plate, but with a known, fixed velocity. The data is all stored in a gigantic text file though, so just extracting it is another project.

In the mean time I'm working on getting the adaptive controller working in simulink with the actual experiment. The last week or so I've spend looking at different ways to put disturbances generated from a SS model on the DM. The most obvious way is to project the desired phase onto some kind of phase poke matrix (maps the actuator inputs to phase, not slopes). In one approach I used the theoretical poke matrix from the manufacturer that is basically a model of the DM surface on a high-resolution grid of points. I don't actually have this for DM61, so in that case I used my estimated poke matrix multiplied by the phase reconstructor. Either one produces phases that seem to flow somewhat like the original model, although I haven't done any kind of analysis to verity that.

As an example of my precise control of Matlab's video functionality, here's the desired phase profile generated from the state space model, positioned on the same size grid as the WFS measurements



And here's the resulting measurements after doing a least-squares fit and applying the corresponding DM commands



They're indistinguishable I know. The fact that the actual measurements show anything that could be believably described as a "flow" is a major success in my book.

Moving blows.

5.30.10

I've been sidetracked the last couple weeks working on this idea of quantifying frozen flow layers. The problem is to identify how many layers are moving in a phase profile, and estimate their velocities. A common idea among some "predictive" AO controllers is to then use this information to generate control commands some number of steps in the future.

Of course, this only works when the velocities are constant and pretty well known, but it seems to be a common approach among certain fields. Based on some comments I've heard, they like it because it incorporates some knowledge about the physics behind the problem. I think deep down some of them just don't trust the completely black box methods that's common in hard core controls applications.

We typically use one of these feared methods to identify a state space model for the turbulence. One question that's bothered us though, is how can we extract the velocity and layer information? Since the controller developed from the state space is optimal, the velocity info has to be embedded in there, but since the states are a product of the ID it isn't clear how.

One approach people in the AO community have tried is to generate a bunch of image correlations from the data, and look for peaks. If the phase is composed of a finite number of layers moving with distinct velocities, the correlations between images separated by enough delay should develop peaks corresponding to each layer. In our case, we have a state space model. And while we could just generate a sequence of data and use these methods, I'd be cooler if we could identify the velocity straight from the system matrices directly.

Coincidentally, I've been reading this book on subspace identification, and it has a good review about calculating the state and output covariance matrices for a state space system. Its very easy to compute the covariance matrices for any number of time steps, so I started to wonder if you could compute the covariance function directly from these matrices. After much, much head banging, it turns out you can.

I don't want to reveal the exact details, but I'll just say even though its not theoretically complicated, its pretty cumbersome in the 2D case, and required many many cups of coffee and a nontrivial amount of cursing to figure out. I still haven't tried it out on real data, but in all the simple 2 layer, integer velocity cases I've developed it works swimmingly, and seems relatively robust to random similarity transformations to the state space.

5.17.10

Now that I can identify simple disturbance models from open-loop data, I'd like to try designing an optimal controller using the identified system. I think it should work out in simulation, but ideally I'll be able to implement it in the experiment and see some results. This is actually much more simple than using actual wavefronts and doing a multichannel problem; since the plant is particularly simple (with a slow enough sampling time), the optimal controller should do a pretty good job.

1. Make sure models can be id'ed using the disturbance model, probably using the 61 actuator modes constructed using the poke matrix.

2. Ignore the PI controller for now? Try to apply disturbances while the integrator is running and see what the results look like.

3. Come up with a script that calculates the optimal controller structure using the controller SS model. First by solving a WH problem, then also by solving a finite time LQR problem. Results should be the same with each.

4. Apply the controller in simulation and in the experiment if all goes well.

5. ....?

6. Profit.

The one wrinkle in all this it that controller and disturbance DM are using different sets of modes, so I have to think about the best way to unify them. Ultimately, I suspect I'll end up using the DM61 modes for generating disturbances, but since that process is basically opaque from the controller's perspective, all the ID and control will be done on the basis of the DM31 modes.