Filtering with Discrete transfer functions in matlab

[gandrin]

Got a DSP problem. I think this is a bit complicated, and may need some DSP gurus to answer it. I've been banging my head on this literally for months now. Thought I had it figured out, but today find I'm still not done.

I have some Impedance spectroscopy data from electrodes, sampled at discrete frequencies. I want to model the transfer function that these electrodes form when connected in series to an amplifier (a simple R||C cct).

The goal is to multiply that discrete transfer function with the FFT(signal), then take IFFT(Product) to show what the amplifier actually sees.

My first inclination was to do the following:

(note: this is a voltage divider across one impedance [ R_amp||C_amp ] that is in series with Z_electrode)

Z_amp= R/sC / (R + 1/sC)
Z_amp= R / (sRC + 1)
Z_divider=Z_amp/(Z_amp + Z_electrode)
Z_divider = R./(R + (Z_electrode .* (1+sRC)))
so now I took the leap of faith...
w=1:1:100000; % the frequencies at which Z_electrode was sampled
Z_divider = R./ (R + (Z_electrode(w) .* (1 + iwRC)));

Then
Y=fft(signal).*Z_divider
filtered_signal=real(ifft(Y));

This gives a reasonable answer...but I'm worried it's false. The reason why is because when I do the same process with just an RC circuit (in this case a generic low-pass RC config), the answer is not right...

RC=1e-5;
Z_lowpass = 1./ (1 + iwRC);

The resulting "filter" was junk, not holding a DC voltage and "anticipating" voltage changes.

signal= ones(1,100000);
signal(50000:70000)=100;
Y=fft(signal);
ZY=Z_lowpass.*Y;
plot(real(ifft(ZY)))
hold
plot(signal,'r')

So I have two questions.

1. I have figured out the coefficients for that lowpass filter
k=deltaT/ (RC + deltaT)
filter([0 k], [1 -(1-k)], signal)

That was NOT EASY to figure out, though wikipedia gave guidance for this particular cct.

So is there an easier way to do that with other circuits in Matlab?


2. So how do I implement this filtering function in the Z_divider example above? How can I "discretize" my RC cct to be able to voltage divide with my electrode Impedance Spectroscopy? ( I do have an equivalent cct for the electrode, but it's a Constant Phase Element, and essentially impossible to figure out analytically.)

 

[gandrin]

I think I've solved my own problem. The problem was the bizarre way MATLAB does FFT. I redid it with the following:

L=100000;
Fs=###;
w_filter=[-L/2:L/2-1]*Fs/L; % have to evaluate the function from -Fs/2:Fs/2-1

Z_lowpass = 1./ (1 + iw_filter*RC);

signal= ones(1,100000);
signal(50000:70000)=100;
Y=fft(signal);

%here's the key--fix the strange Matlab orientation
Y_shift=fftshift(Y);

ZY=Z_lowpass.*Y_shift;
ZY_unshift=ifftshift(ZY);
plot(real(ifft(ZY_unshift)))
hold
plot(signal,'r')

 

[JJacquelin]

I do have an equivalent cct for the electrode, but it's a Constant Phase Element, and essentially impossible to figure out analytically.

Hi !
On a theoretical viewpoint, CPE can be figured out as shown in the paper "The Phasance Concept" published on Scribd.
For practical applications, many equivalent networks can be derived and used to model real electrical components. For example: "Theoretical impédances of capacitive électrodes" published on Scribd :
http://www.scribd.com/JJacquelin/documents

 

[rbj]

the thing you got to remember about the DFT (this is not anything about MATLAB per se) is that it has an inherent periodicity. so the last half of the data (which precedes the first half of the repeated cycle of the same data) can be thought of as the "negative frequency" or the "negative time" data that precedes the first half. MATLAB's fftshift() is just meant to swap the last half with the first so that it looks like the negative frequencies precede the positive frequencies.the other dumb thing MATLAB does is insist that all indices are positive integers. 0 or negative integer indices are not allowed in MATLAB which is unfortunate since there is so much literature, especially in DSP, that has non-positive indices.