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.
Best Answer
From your description of the experiment (please correct me if my assumptions are wrong), it sounds like your apparatus consists of the application of a controlled stress to the sample (and the sensor), and the resulting strain in the sensor is measured. Whenever the stress applied by your apparatus changes, it will take some time for the system to settle to it's new equilibrium. It could be that sampling at $1 Hz$ is allowing plenty of time for equilibration, but sampling at higher frequencies you are recording the oscillations of the system as it has not yet settled.
One way to test would be to run the experiment without changing the applied force, just recording the strain at various sampling rates, and looking to see if the noise spectrum still depends on the sample rate in the way you describe. If it does, then the noise is a result of the frequency dependence of the electronics. If it does not, then the noise is resulting from the physical behavior of the sample