White Noise Generator Through Low Pass Filter - Autocorrelation of normalized data dropping below 0 (showing negative correlation)

Question

I am doing the following experiment: http://physlab.lums.edu.pk/images/a/ab/Correlation.pdf

It is basically a white noise generator running through a low pass filter. Sample voltages are taken and normalized, then autocorrelated using a few methods. The autocorrelation using normalized data is supposed to start at a value of 1 at 0 lag (which it does) and then exponentially decay to 0(which it doesnt). However it is decaying below 0.

Also The auto generated coefficient for an exponential fit (a*np.exp(-(t/b)) + c) (b should = RC), does not match the RC constant of the Low pass filter (R=10kohm, C=0.5 microF, RC=5ms).

The methods I am using for the autocorrelation are as follows: f=list of voltage values with mean subtracted from each value, 100,000 samples @ 10000 S/s

Using Built in Correlate function:

def autocorr(f):

    result = np.correlate(f, f, mode='full')

    return result[result.size/2:]

autocorr=autocorr(f)

autocorr /= autocorr[autocorr.argmax()]

Using FastFourierTransform:

N = len(f)

fvi = np.fft.fft(f) #n=2*N)

acf = np.real( np.fft.ifft( fvi * np.conjugate(fvi) )[:N] )

Plotting using inbuilt python autocor plot:

import matplotlib.pyplot as plt

plt.acorr(f)

Which you can see does not drop below 0. I cannot find out the exact code that is executed with the plt.acorr(f) as it does not support outputting of data.

Here is the fitted exponential (clearly showing dip below 0) fitting function =

from scipy.optimize import curve_fit

def func(t,a,b,c):

return a*np.exp(-(t/b)) + c

popt, pcov = curve_fit(func, t, acf)

However the autogenerated constant from the exponential fit, that should equal RC (5ms) actually comes out as 1.043822ms Also when I reduce the RC by 50% the calculated coefficient reduces by 33%. Also if I change the sampling frequency this reduces the coefficient

Is there a way to produce the autocorrelation of the data gathered, that consists of an exponential decay to 0 with the decay constant as the RC of the low pass filter circuit?

I have tried various RC values for the LPF, they all dip below 0, in fact the lower the RC, the lower it dips below 0.

Could this be related to an incorrect offset voltage of the white noise generator? I changed the offset, and it did increase the lowest value that the autocorrelation drops to, but I am not sure if I adjust it again it will fix the problem completely.

So to recap, 2 issues

• Autocorrelation does not decay to 0, it drops below 0

• Fitted exponential curve decay coefficient, does not equal the RC constant of the circuit.

This thread Autocorrelation of noise - negative correlation Seems to be related to the same experiment and the same issue (unanswered)

                      **UPDATE 7-3-15:**

Things that have been tried: The White Noise generator has a threshold knob to set the threshold frequency. It did not fix either of the issues: Highest Threshold Setting:

Lowest Threshold Setting:

Autocorrelation was fitted to Python generated White noise filtered through simulated RC low pass filter circuits (different RC values). This produced correct results. This would lead me to believe the code being used to autocorrelate is correct.

I have tested the FFT method and inbuilt python correlate function to autocorrelate data and they give the same results, also confirming the code should be ok.

A physical signal generator, generating square waves was tested through the A2D and plotted v/s time. The result was correct with the correct timings matching the sampling frequency. This confirms the A2D is functioning correctly.

Example of square wave through A2D:

The square waves were also ran through 2 RC low pass filter circuits. The correct RC values were calculated using an exponential fit. This confirms the RC circuits are functioning correctly and also the code used to fit the exponentials is functioning. However decay below 0 occurred

Example of fitted exponential on filtered square wave:

There is a voltage offset in the white noise generator, this was adjusted a few times (above and below 0V), it did not fix the incorrect fitted exponential RC constants. It does effect the decay of the autocorrelated data, but not to the point where it can be determined if this is the cause of the issue.

Plot below. B=before 1st offset change, A = after 1st offset change, A2-x = after 2nd offset change. RC's correspond to eachother e.g. B1, A1, A2-1 have the same RC etc

I have checked the power spectra of 2 samples of non filtered noise from the white noise generator with the following results:

Please ignore the units

It is possible the noise being generated isn't white noise at all. I will need to check with more samples.

The data does fit autocorrelated white noise (decaying exponential), the problem is, it decays beyond 0 and then recovers and the decay constants are roughly a factor of 10 out from the actual RC values of the low pass filters that the white noise is being filtered through. Currently testing the power spectra of the white noise (not filtered or normalized).All plots of data are just samples with a few samples, I have ran more, to many to include. Thanks to User:Floris for suggestions on what to look at.

                      **UPDATE 7-15-15:**

Well it turns out the problem does lie with the white noise generator. A very strange issue.

I added a voltage follwer circuit before the RC circuit (to increase the current available) and the results improved dramatically. I'm not sure if this is the final fix as there is some weird issues with the voltage follower, adding voltage gain for some reason. I tested with another white noise source before using the voltage follower and I was able to achieve ~ 90% accuracy comparing the exponentially fit RC value to the theoretical RC value and can now match that with the original white noise generator and voltage follower circuit.

Example of a collection of 10 sample sets @ RC=5ms:

Compared to the "working" white noise source:

I will update with a final answer when I know exactly whats happening and why the issue was occuring. If anyone has any suggestions also, much appriciated.

Answer 1

A single measurement like this has a lot of noise on it - and random signal is always going to have some random correlation. You should definitely not pay too much attention to the stuff that is in the tail of the correlation distribution - it's all noise.

The fact that the built in function does not produce negative values is related to you only looking at the first 10 bins (which is the default). If you plotted the full range, it would go negative just like the others - use

plt.acorr(f, maxlags=None)

Incidentally, it does return all the values in four variables - if you do

result = plt.acorr(f, maxlags=None)

you can examine result[0] and result[1] to give you the values the function calculated - you will find it's similar (identical) to the one you get using other functions. The reason the plot is different, as I said, is that you are not looking at the same range of lags.

Let me demonstrate a few principles. I will be generating some white noise, filtering it with the response of an RC circuit, and plotting the autocorrelation multiple times. We are going to learn how to create white noise (and how to filter it); that the autocorrelation will not be the same twice; that it can go negative, even if you take the mean of 20; and that the measured time constant will vary by quite a bit. Bottom line - you will come away with a better feeling for this. There is much more to learn...

Here is the code:

# white noise

import numpy as np
import matplotlib.pyplot as plt
import math
from scipy import stats

def autocorr(f):
    temp = np.correlate(f, f, mode='full')
    mid = temp.size/2
    return temp[mid:]

plt.close('all')

# simulate a simple RC circuit
# V = 1/jwC/(R+1/jwC) = 1/(jwRC+1)
C = 1e-6      # 1 uF
R = 1.0e3     # 1 kOhm
RC = R*C
fs = 100000   # sampling frequency
fc = 1.0 / RC # cutoff frequency
repeats = 100 # number of times experiment is repeated
ns = 10000    # number of samples per repeat (1/10th of a second)
nplot = 500   # number of samples in autocorrelation to plot

t = np.arange(ns)/(1.0*fs) # time corresponding to each sample
fr = fs / (2.0 * ns) # resolution per bin in the FFT

freq =  np.arange(ns) * fr  # frequency bins in FFT
response = np.divide(1, 1j * 2 * math.pi * freq * RC + 1) # complex response per bin

# plot response function:

fn = 4* np.floor(fc / freq[1]) # two octaves above 3dB point
plt.figure()
plt.subplot(2,2,1)
plt.plot(freq[:fn], np.abs(response[:fn]))
plt.xlabel('frequency')
plt.title('Filter response: amplitude')

plt.subplot(2,2,2)
plt.plot(freq[:fn], np.angle(response[:fn], deg=1))
plt.title('Filter response: phase')
plt.ylabel('angle (deg)')
plt.xlabel('frequency')

fa = []  # array where we will store the filtered results
plt.subplot(2,2,3)
for j in range(repeats):
    e = np.random.normal(size=ns)
    # take Fourier transform:
    ff = np.fft.fft(e)
    # roll off the frequencies
    fr = ff * response
    # take inverse FFT:
    filtered = np.fft.ifft(fr)
    ac = autocorr(filtered)
    acNorm = np.abs(ac / ac.max()) # <<< using abs function here: get amplitude
    fa.append(acNorm)
    plt.plot(acNorm[:nplot])

plt.xlabel('correlation distance')
plt.title('Autocorrelation: result of %d samples'%repeats) 
plt.show()

# calculate the mean of the repeated samples:
repeatedMean = np.mean(fa, axis=0)

plt.subplot(2,2,4)
plt.semilogy(repeatedMean[:nplot])
plt.title('mean of %d repeats'%repeats)
plt.xlabel('correlation distance')
plt.show()

# fit an exponential to the data
# anything beyond the first 20 points does not make much sense - contains no useful information

# fit an exponential to the data
# anything beyond the first 20 points does not make much sense - contains no useful information

nfit = fs*RC  # <<<< added
slopes = []
plt.figure()
for j in range(repeats):
    slope, intercept, r_value, p_value, std_err = stats.linregress(t[:nfit], np.log(fa[j][:nfit]))
    plt.plot(t[:nfit], fa[j][:nfit])
    slopes.append(slope)
plt.title('autocorrelation used for fitting') # <<< fixed typo
plt.xlabel('correlation time')
plt.show()
print 'mean time constant is %.2f ms'%np.mean(np.divide(-1000,slopes)) # <<< changed
print 'standard deviation is %.3f ms'%np.std(np.divide(1000, slopes))  # <<< changed

And here are the plots it produces:

The time constant derived from these fits (inverse of slope) averaged over 100 fits was 1.01 ± 0.15 ms - consistent with the time constant of 1.0 ms that was implied in the choice of components in my code (1 kOhm and 1 uF).

Play around with the code - and let me know if this clears things up for you.