Low Pass Filtering - Generate a Slowly Changing Signal in Matlab

[sarahpaninsky]

Hi,
I want to generate a slowly changing signal (1 second duruation).
For this, I must filter white noise with a filter which is characterized by an alpha function as below

alpha(t)= t*exp(-t/tau) ( I think function is low-pass filter response)

where tau=3 ms. I couldn't achive this in matlab. Can anyone help me for this?
Regards

 

[marcusl]

Can you show us your work? We can help you better if we see your approach and your Matlab code.

By the way, is your alpha a step response (probably) or an impulse response?

 

[sarahpaninsky]

Can you show us your work? We can help you better if we see your approach and your Matlab code.

By the way, is your alpha a step response (probably) or an impulse response?

clear
figure
tau=3;
noise=8+7*randn(1,1e5);
t=[0:0.01:1e3];
y=conv(noise,t.*exp(-t./tau));
subplot(3,1,1)
title('Gaussian Noise')
plot(noise)
subplot(3,1,2)
plot(y)
title('Filtered Noise')

My code is like above.
As I said before the filtered signal must be slowly changing, is like a stimulus not actually like noise to a biological neuron. Please modify my code

 

[marcusl]

Be careful of your time units. You specified tau=3, implying that your waveforms are sampled at 1000 Hz. Then each point in your sequence is separated by 1 ms. Your impulse response dies to essentially zero after 30 ms, so there's no need to include 1e5 points!

Below I relabeled your t as n to make the point that this is the index of a sampled data sequence, and reduced the length of the impulse response h to 30 points. I also reduced the input sequence noise to 10,000 points so it runs faster, but this is not important.

clear
tau=3;
noise=randn(1,1e4);
n=[0:30];
h = n.*exp(-n/tau);
y=conv(noise,h);
figure
subplot(3,1,1)
title('Gaussian Noise')
plot(noise)
subplot(3,1,2)
plot(y)
title('Filtered Noise')

This does low-pass filter your noise. You can display the frequency response of your filter and output with the following commands:

Fs = 1000;
Hs=spectrum.periodogram
figure
psd(Hs,h,'Fs',Fs)
title('Filter Response')
figure
psd(Hs,y,'Fs',Fs)
title('Output Spectrum')

Fs above is the sampling frequency, 1000 Hz in this case. Notice that your cutoff frequency (-3 dB) is about 35 Hz. You probably studied the relation between time constant and cutoff frequency fc; what is it? See if you can adjust tau to get fc=1 Hz. You'll have to readjust the sampling frequency to match your new tau...

 

[sarahpaninsky]

Be careful of your time units. You specified tau=3, implying that your waveforms are sampled at 1000 Hz. Then each point in your sequence is separated by 1 ms. Your impulse response dies to essentially zero after 30 ms, so there's no need to include 1e5 points!

Below I relabeled your t as n to make the point that this is the index of a sampled data sequence, and reduced the length of the impulse response h to 30 points. I also reduced the input sequence noise to 10,000 points so it runs faster, but this is not important.

clear
tau=3;
noise=randn(1,1e4);
n=[0:30];
h = n.*exp(-n/tau);
y=conv(noise,h);
figure
subplot(3,1,1)
title('Gaussian Noise')
plot(noise)
subplot(3,1,2)
plot(y)
title('Filtered Noise')

This does low-pass filter your noise. You can display the frequency response of your filter and output with the following commands:

Fs = 1000;
Hs=spectrum.periodogram
figure
psd(Hs,h,'Fs',Fs)
title('Filter Response')
figure
psd(Hs,y,'Fs',Fs)
title('Output Spectrum')

Fs above is the sampling frequency, 1000 Hz in this case. Notice that your cutoff frequency (-3 dB) is about 35 Hz. You probably studied the relation between time constant and cutoff frequency fc; what is it? See if you can adjust tau to get fc=1 Hz. You'll have to readjust the sampling frequency to match your new tau...

Thanks for your interest Marcusl,
But the results are not sufficent.
I want to show you my work, I want to reproduce the low-pass gaussian white noise stimuli as in the image I attached. if you look at it, you may understand my problem with the code. What could I do ? Filtering is Okey but the the code is fails for amplitude and fluctuating?

 

[marcusl]

Sorry, I thought this was a homework assignment. I'll take a close look at this tonight.
What are the spikes they refer to? Are they individual channel responses? What is the significance of jitter? (I work in radar, not membranes...)

Would you post the paper, or link to it?

thx

 

[sarahpaninsky]

my e-mail is "opamp67@hotmail.com", if you send a blank mail, I will send you immediately. Because I cannot attach the paper as the file size limit of the forum.

This is an scientifistic paper on Neuron modelling in biomedical disciplinary. in the paper Figure1 is important for me, especially top two figures which are related to low-pass gaussian noise. My final code is below to produce the signal Fgi1B (top trace)


clear
figure
tau=3;
noise=8+7*randn(1,5e4);
t=[0:0.01:499.99];
h=t.*exp(-t/tau);
y1=(fft(noise));
y2=(fft(h));
y3=(ifft(y1.*y2));
plot(t,y3);
axis([0 500 -2 11])

 

[marcusl]

I can take a few moments at lunch to make the following comments:
1) You've used tau = 3 sec in your code instead of 3 msec. This will produce a very low frequency output (filter cutoff of 0.05 Hz, according to calc below).

2) You can predict the cutoff (-3dB or half power) frequency of the analog alpha filter as follows. You've specified an impulse response
$$\alpha(t)=t\exp{(-t/\tau)}$$ which has Laplace transform
$$F(s)=\frac{1}{(s+1/\tau)^2}$$.
Evaluating $$s=i\omega+\sigma$$ on the frequency (imaginary) axis
gives the frequency response as
$$|F(\omega)|=\frac{1}{|\omega+1/\tau|^2}$$.
This is a two pole lowpass filter with cutoff frequency
$$f_c=\frac{1}{2\pi\tau}$$.
The point is, $$\tau=3 ms$$ gives a filter rolloff of $$f_c=53 Hz$$.This is roughly consistent with the waveform in Fig. 1B, and is not consistent with your statement that the waveform should change on the order of seconds.

3) You can see from the expression for F that the filter gain is not normalized to one at zero frequency. Dividing by tau^2 gives proper DC normalization
$$\alpha(t)=\frac{t}{\tau^2}\exp{(-t/\tau)}$$.

4) In general, sampling the impulse response of an analog filter the way you've done it here is not the optimal way to implement a digital filter. One problem is that analog and digital cutoff frequencies don't match. If you look at the power spectral density of h using the psd commands I gave above, you'll see that the -3dB frequency response of the digital filter using tau=3 ms is around 35Hz instead of the 53 Hz calculated for the analog filter. A sampled impulse response implementation also aliases, so F(w) doesn't roll off towards zero as expected but flattens out. Finally, the DC normalization is thrown off. I assume that is why Fig. 1 shows a mean of about 4, even though the input had a mean of 8 that should be unchanged by the lowpass.

I'm not sure whether these considerations are important to your paper and work, but at least you have an appreciation of what this system is doing from the standpoint of digital signal processing.

 

[marcusl]

Sorry, there's a typo in the frequency response--there should be an i in front of omega:

$$|F(\omega)|=\frac{1}{|i\omega+1/\tau|^2}$$

Of course this simplifies to
$$|F(\omega)|=\frac{1}{\omega^2+1/\tau^2}$$

 

[marcusl]

Ok, the "scientifistic" paper you sent contains no information on the signal processing used within, beyond that mentioned in the Figure attachment above. The fact that the signal mean values in Fig. 1 A,B differ from input to output is strange. We'll take the "filter method" in Fig. 1 at face value, however, lacking any other information. Changing the OP's code as follows at least approximates their results:

clear
tau=3e-3;
noise=8+7*randn(1,1e4);
t=[0:0.001:0.03];
h = t .* exp(-t./tau) ;
y=conv(noise,h/sum(h));
figure
subplot(3,1,1)
plot(noise)
title('Gaussian Noise')
axis([500 1000 0 15])
subplot(3,1,2)
plot(y)
title('Filtered Noise')
axis([500 1000 0 15])
figure
plot(h)
title('Filter Impulse Response')

One-half second of signal is plotted, as in the paper. The filter is normalized ad hoc by dividing by the sum of the impulse response.

 

[sarahpaninsky]

Ok, the "scientifistic" paper you sent contains no information on the signal processing used within, beyond that mentioned in the Figure attachment above. The fact that the signal mean values in Fig. 1 A,B differ from input to output is strange. We'll take the "filter method" in Fig. 1 at face value, however, lacking any other information. Changing the OP's code as follows at least approximates their results:

clear
tau=3e-3;
noise=8+7*randn(1,1e4);
t=[0:0.001:0.03];
h = t .* exp(-t./tau) ;
y=conv(noise,h/sum(h));
figure
subplot(3,1,1)
plot(noise)
title('Gaussian Noise')
axis([500 1000 0 15])
subplot(3,1,2)
plot(y)
title('Filtered Noise')
axis([500 1000 0 15])
figure
plot(h)
title('Filter Impulse Response')

One-half second of signal is plotted, as in the paper. The filter is normalized ad hoc by dividing by the sum of the impulse response.

That is GREAT Marcusl. I am glad to you. Everything is OK now. Thanks your kind interest.