DSPRelated.com
Forums

Inverse filter

Started by Rahul02111997 6 years ago9 replieslatest reply 6 years ago287 views

Hello Everyone,

I am passing a signal through a RC HPF (cutoff freq = 14000 Hz). The sampling rate at which the filtered signal is acquired = 400 kS/s. 

Now suppose, I don't know my input signal and only have the output signal from the output terminal of this HPF. How can I get the time domain signal which must be fed into the filter to obtain the given output signal? 

For getting the PSD of the input signal,  I am simply dividing the periodograms by the square of magnitude response. It works perfectly. If my input is white noise into the HPF. After dividing the Periodogram of the filtered signal, I am easily getting the periodogram of white noise i.e. flat spectrum. But I need to get the time series data which must be input to the system. How can I do this? 

In short, I want to reverse the effect of the HPF to get the input time domain data. 

Thanks in advance!

[ - ]
Reply by Rick LyonsJune 13, 2018

Hi Rahul. I think you want to perform what's called "deconvolution." And deconvolution can be performed (I) using your filter's output sequence (time-domain deconvolution) or (II) using the spectrum of your filter's output sequence (frequency-domain deconvolution).

If you have MATLAB, type "help deconv" in the Command window and then click you mouse on blue "Reference page for deconv" letters. Then you'll see the following example of time-domain deconvolution:

deconvolution_12969.jpg

In the above example, translated to your signals, we have:

  u = the time-domain impulse response of your filter (You have this.)

  v = the filter's input sequence (WHAT YOU WANT!)

  c = the output sequence of your filter. (You have this.)

Rahul, try out the above time-domain 'deconv()' command and let us know if it works for you. If you would like to discuss the notion of frequency-domain deconvolution, let us know.

[ - ]
Reply by Rahul02111997June 13, 2018

Thanks a lot Rick!

It works perfectly and is quite intuitive way. I was not aware of the this deconv function. Thanks a lot for your suggestion. Can you also comment on the stability of this approach? Also, can you tell me why inverting the poles and zeros of the filter technique do not work as good as this technique? I guess in both methods, we are doing similar processing, we are just removing (deconvolving) the filter response from the filtered data

Just for reference to others who might have same doubt in future, here is the simple Matlab simulation code to verify this method: 

clear all

close all

fs = 5000;

L = 5000;

R = 1000;

t = 0:1/fs:60-1/fs;

data = randn(length(t),1);

[pxx,freq] = pwelch(data,L,R,[],fs);

figure(1)

plot(freq,10*log10(pxx), '--');



%% filtering using butterworth filter

order    =1;

fcutlow  = 1000;

fcuthigh = 15915;

% [b,a]    = butter(order,[fcutlow,fcuthigh]/(500000/2), 'bandpass');

[b,a] = butter(order,155/(fs/2),'high');

sig0_filt  = filter(b,a,data);

[pxxf,freq1] = pwelch(sig0_filt,L,R,[],fs);

figure(2)

plot(freq1,10*log10(pxxf), '--');



%% filtering in time domain using convolution

[hh,tt] = impz(b,a);

Time_filtered = conv(data,hh);

[pxxTDHPF,freq1] = pwelch(Time_filtered,L,R,[],fs);

figure(3)

plot(freq1,10*log10(pxxTDHPF), '--');



%% Reconstruction of unfiltered data by deconvolving with filter impulse response

recons_dataa = deconv(Time_filtered,hh);

[pxxTDrecons,freq1] = pwelch(recons_dataa,L,R,[],fs);

figure(4)

plot(freq1,10*log10(pxxTDrecons), '--');

[ - ]
Reply by Rick LyonsJune 13, 2018

Hi Rahul. Sorry for my delayed reply. I haven't worked on any projects that required me to study deconvolution in any detail. After searching the web a little I have the distinct impression that deconvolution of noisy signals may be a bit complicated.

Also, I don't know what you meant when you wrote: "...inverting the poles and zeros of the filter technique..."

[ - ]
Reply by achesirJune 13, 2018

Design a filter whose transfer function is the reciprocal of the filter you have. Feed the output of your current filter to the input of this second filter. The output of the second filter will be a copy of your input signal.

[ - ]
Reply by kazJune 13, 2018

you mean get LPF to filter a signal that has only power left after HPF.

Are you sure as I imagine nothing is left? and How can we recover filtered section

[ - ]
Reply by achesirJune 13, 2018

Rahul stated that he does not know his input signal, but he presumably knows the design of the HPF. Therefore it should be straightforward to determine the transfer function of his HPF. Therefore, a filter with a reciprocal for the HPF, which will require some active component to apply the reciprocal gain in the stop bands of the HPF, should be realizable, albeit with some delay.

[ - ]
Reply by kazJune 13, 2018

well with fixed point data that is impossible in my opinion. But may be it is true for floating point.

[ - ]
Reply by jbrowerJune 13, 2018

Kaz-

I've calculated:

  H(s) = Y(s)/X(s)

many times before before but to obtain reliable results I usually had to do two things (i) use 32-bit fixed-point FFT results, and (ii) average f-domain frames, maybe 3-6 or so.

Performing a similar calculation for X(s) and then finding x[n] should be a do-able thing, once you have come up with a reliable, consistent method to calculate X(s) based on your system conditions / operating parameters.  As Rahul mentions, to obtain continuous, phase-coherent x[n] across f-domain frames there are other issues, but I think you are rightly worried about precision and other practical aspects, so first focusing on X(s) seems like a good idea.

-Jeff

[ - ]
Reply by Rahul02111997June 13, 2018

Thanks for the reply!

Can you give me more detail about this? May be a small Matlab code or some link? 

Because I guess it is not that straightforward. There will be problem at sampling frequency/2 as stated in this forum:

https://dsp.stackexchange.com/questions/24557/how-...