How to filter a signal with lowpass/highpass?

Hi all,

I am completely new in this field. Therefore, I need some help from the scratch.

According to a protocol, I have to filter a signal at 1.2 kHz ("the current derivative (of the signal) was addicionally filter at 1.2 kHz (sic)).

I've tried to filter the signal using the funtion lowpass:

 y = lowpass(x,fpass,fs) %specifies that x has been sampled at a rate of fs hertz. fpass is the passband frequency of the filter in hertz)

x is my data, fpass is 1200 (1.2 kHz), but fs? I've checked the protocol and I've found this: "current signals were filtered at 2 kHz and digitized gap-free at 25 kHz."

So, what is the fs in my case? 2000, 25000, or 27000 (2000 + 25000)?

My main question is: how can I filter the signal at 1.2 kHz?

I've attached the data (Trial file). There are two variables (1x26 double), Time (s) and Intensity (pA). Plot it as (Time,Intensity).

Many thanks,

Jose

Answer 1

[Star Strider]

The easiest way to determine the sampling frequency is to measure it:

D = load('trial1.mat');
A = D.Intensity;
t = D.Time;

Ts = mean(diff(t));                             % Sampling Interval
Fs = 1/Ts;                                      % Sampling Frequency

The sampling frequency is 6250 (Hz, assuming ‘t’ is in seconds), and appears to have been sampled regularly (i.e. the sampling intervals are essentially identical). Define the passband of a filter so long as it is less than the Nyquist frequency, Fs/2, so a passband or stopband of 1.2kHz is possible.

I have no idea what ‘...current signals were filtered at 2 kHz and digitized gap-free at 25 kHz.’ means. You have to sort that.

[Jose Rego Terol]

Sorry Star Strider. I cannot give you more than "one vote" and "accept this answer".

This code is a piece of cake. Your explanation was straightforward even for me.

Many, many thanks! It works!

Now, how can I learn what you did here? Where can I start?

[Star Strider]

As always, my pleasure!

I first use the load function to bring the contents of ‘trial1.mat’ into the workspace. Hovering over it, the Editor tooltip told me that it had two fields, so I assigned each of those fields to a variable.

The diff function takes the successive differences of ‘t’, so taking the mean of those defines the sampling interval, ‘Ts’. (I also calculated the standard deviation of diff(t), that was on the order of system floating point precision, or essentially zero, so the sampling interval is regular.) The sampling interval is in seconds (assuming that ‘t’ is in seconds), so the sampling frequency ‘Fs’ is the inverse of that. The Nyquist frequency is the highest frequency in a sampled signal that can be uniquely identified.

If you have any other questions, I will do my best to provide answers.

[Jose Rego Terol]

"The diff function takes the successive differences of ‘t’, so taking the mean of those defines the sampling interval". I didn't know that. Thanks for the explanation :)

[Star Strider]

As always, my pleasure!

It actually doesn’t define the sampling interval. (That was an unfortunate misstatement on my part.) The sampling intervals are defined by the experiment design and the ADC hardware. The diff function detects the sampling intervals and permits estimation of the sampling intervals regularity (using the std and related statistical functions).

[Jose Rego Terol]

By the way, as the protocol says "the current derivative (of the signal) " I was using diff for this purpose.

z=diff(I). But then I realised that I need to include the time in this equation (doc diff) to get the first derivative of the current (see section Approximate derivatives with diff). For that reason, I changed the equation as follows: The recording frequency is 0.00625 micro Hz (one point every 160 microseconds, 1/160). Therefore the diff formula is: z=diff(I)/0.00016. Although the signal pattern looks fine, the magnitude of the data is too high.

Therefore, I googled for the solution, and I have found this function, gradient. I have applied this function to the Intensity (one of the fields of the file). Both gradient (I) and diff(I) looks pretty similar. The point is, in what should I rely on to get the first derivative of a signal?

[Star Strider]

‘The point is, in what should I rely on to get the first derivative of a signal?’

Definitely use the gradient function. The principal reason is that it is designed to calculate the numerical derivative, so the output is the same size as the input (as opposed to the diff function that simply takes the successive differences in the vector, and that reduces the size in the chosen dimension by at least 1).

The gradient function wants a constant sampling interval, that being the second ‘h’ argument. However it is possible to get around that by doing element-wise division of the gradient of the function and the gradient of the independent variable vector (here: time):

dydt = gradient(y) ./ gradient(t);

giving the time derivative corresponding to the sampling times.

I am not certain that I understand your sampling interval and frequency derivations.

With:

Ts = 160E-6                                 % Sampling Interval (sec)
Fs = 1/Ts                                   % Sampling Frequency (Hz)
Fn = Fs/2                                   % Nyquist Frequency (Hz)

the appropriate results (as I calculate them) are:

Ts =
                   0.00016

Fs =
                      6250

Fn =
                      3125

[Jose Rego Terol]

It is crystal clear that function gradient fits better in my situation. Thanks for the confirmation.

Regarding the sampling interval, I calculated it in this way

t(2)-t(1)

Looks naive, but the result is the same than yours. Hence, the rest of the parameters are also the same. Your way is fancier than mine, though.

About gradient function input: The gradient function wants a constant sampling interval, that being the second ‘h’ argument. I expected 'h' argument as a number, not a vector, since MatLab doc for gradient uses the sampling interval as 'h' argument. You claim that the element-wise division of the gradient of the function and the gradient of the independent variable vector gives you the 'h' argument. However, dydt is a column vector. Hence, I cannot uses dydt as an argument for gradient function, right?

Another question that arises is that I do not understand for what purpose you are applying this equation

dydt = gradient(y) ./ gradient(t);

having the gradient function

z = gradient(y,Ts)

I've calculated the mean and the SD of both z and dydt. They are pretty the same:

mean(dydt)=

       -0.0767719681812023

mean(z)=

       -0.0767719680460206

std(dydt)=

          13531.3409619448

std(z)=

          13531.3409619518

Is there something that I am misunderstanding?

[Star Strider]

The ‘h’ argument is a constant scalar.

‘You claim that the element-wise division of the gradient of the function and the gradient of the independent variable vector gives you the 'h' argument.’

I hope I did not give you that impression. That gives you a reasonably accurate derivative in the event that the sampling intervals differ significantly enough so that one ‘h’ argument would not be accurate. (The default value for ‘h’ is 1, and that works in this instance.)

‘However, dydt is a column vector. Hence, I cannot uses dydt as an argument for gradient function, right?’

For vectors there should be no problem. For matrices, the orders matter, because gradient then calculates the spatial gradient, not simply the numerical derivative, so there are as many matrix outputs as there are matrix argument dimensions.

These are the same, except that they are tranposes of each other:

V = rand(10,1);
dV = gradient(V);
dVT = gradient(V.');

‘Is there something that I am misunderstanding?’

I doubt it. This instance has a regularly-sampled vector, so you could use one scalar ‘h’ without problems. I was demonstrating a situation where the sampling intervals are not regular, so doing the element-wise division produces an accurate numerical derivative that would not be possible otherwise. That does not apply to your vector, I was simply iillustrating that solution.

[Jose Rego Terol]

I doubt it. This instance has a regularly-sampled vector, so you could use one scalar ‘h’ without problems. I was demonstrating a situation where the sampling intervals are not regular, so doing the element-wise division produces an accurate numerical derivative that would not be possible otherwise. That does not apply to your vector, I was simply iillustrating that solution.

Got it!

I want to ask you for a completely different thing. I have a large number of .txt files in a folder and, I want to know which one has a variable (fi spor._M). For example, files a, b, c, and f have these variables but d, and f do not have this variable.

I was using find(fi spor._M) but I get the number of matches (e.g. 8). I do not know how many times this variable is written in each file. It could be that the variable appears in eight files or 1 file, but eight times. Do you catch me?

What I want is to display the name of the files that this variable appears.

[Star Strider]

To the best of my knowledge, it would be necessary to load all the files and then search each file for the variables of interest. I doubt that there is any other way to determine that.