# Reverse engineering a digital filter

I have 2 simultaneous signals that both are designed to measure eye movements. They are sampled at 250 Hz. We have 12 subject recordings. For 3 stable periods in each subject, we choose 256 points and did an FFT. Prior to the FFT, the data were mean-centered and detrended with a 2nd order polynomial. They were also windowed with a Hann window. We are focused on the magnitude spectra plots.
We have a total of 12 X 3 = 36 magnitude spectra, which we average. These averages are shown in the attached figure.

My hypothesis is that the second signal is a low-pass filtered version of the first. The filter has ringing in the passband.

I want to design a digital filter that I would apply to signal 1 that, after fft, analyses would produce a magnitude spectra like that of signal 2.

• The question is which features of the filter are important to you and which aren't. Is the exact behavior in the stop band relevant? I would guess that it is mainly the cut-off frequency and the minimum stop band attenuation that are relevant. – Matt L. Aug 13 '19 at 18:43
• Thanks Matt L. Ideally, everything. But yes, if I had to choose, I would emphasize the cutoff-off-frequency and the minimum stop band. – Lee Friedman Aug 14 '19 at 14:14

Signal 1 was an average of a left and right signal and signal 2 was a binocular signal.

I followed the suggestion of MBaz, and computed the frequency response of the filter for each segment (N = 3) for each subject (N=12). Then I averaged the frequency responses of the filters. Here is the result:

MBaz and Matt: Thank you so much for solving my problem. Matt, it looks close to your filter. Lee

• Now you know the desired magnitude response of the filter, which you can use to actually design the filter. Or try the one in my answer. – Matt L. Aug 16 '19 at 8:28
• Hi Matt, MBaz's approach gives me the filter. Once I have that, I use MATLAB freqz for the frequency response. I think the frequency response is all I need. All I need to say for my manuscript is that it looks like the binocular signal is a filtered version of the LR Average signal and show the frequency response of a filter that would do it. I think I have what I need. What I would really like to know now is: how and why does this filtering happen? We have no schematics and we don't have human readable code. The device is no longer being manufactured. The whole company was bought out. – Lee Friedman Aug 17 '19 at 13:08
• The binocular signal is superior on a number of data quality metrics but this filter removes or obscures some very important fast features of some very small eye movement events. In certain contexts these events are extremely important. – Lee Friedman Aug 17 '19 at 13:15
• OK, it was because in the question you said that you wanted to design a digital filter. – Matt L. Aug 17 '19 at 20:26
• Matt, you are correct. I guess I misspoke. I am completely satisfied with the frequency response for the filter. I will use your words to describe it as a low order low pass filter with marked ringing in the stop band. Am I correct in assuming that the filter could be analog or digital? And is it appropriate to cite the 3db point and the 6db point. – Lee Friedman Aug 18 '19 at 21:12

You want a filter $$h(t)$$ such that $$y(t) = x(t) \ast h(t),$$ where $$x(t)$$ is your signal 1 and $$y(t)$$ is your signal 2.

To find the filter, note that $$Y(f) = X(f)H(f)$$, so $$h(t) = \text{IFFT} \lbrace Y(f) / X(f) \rbrace.$$

• Thank you. This looks like an excellent approach and exactly what I am looking for. Give me a few days to test it out and I will get back to you if I have any problems. – Lee Friedman Aug 15 '19 at 0:13

From looking at the plot, I would assume that the filter is a relatively low order FIR filter designed using the window method or a least squares design criterion.

If you use Matlab/Octave you could try the function firls for a least squares design. I used my own routine lslevin to design the following length $$12$$ linear phase FIR filter:

N = 12;
w = pi*[linspace(0,1/25,50),linspace(3/25,1,1000)];
D = [exp(-j*w(1:50)*(N-1)/2),zeros(1,1000)];
W = [10*ones(1,50),ones(1,1000)];
h = lslevin(N,w,D,W);


I would expect that filtering your first signal with that filter should give a result similar to your second signal. If not, playing around with the parameters should get you relatively close.

• Matt, I sure do appreciate your efforts and I will definitely give it a try. Thank you for your considerable effort. – Lee Friedman Aug 15 '19 at 0:15