CODEWITHSUNDEEP
pythonnumpyiterator
Mish JPFD
Mish JPFD

Reputation: 55

How to calculate Delta F / F using python?

I've recently "taught" myself python in order to analyze data for my experiments. As such I'm pretty clueless on many aspects. I've managed to make my analysis work for certain files but in some cases it breaks down and I imagine it is a result of faulty programming.

Currently I export a file containing 3 numpy arrays. One of these arrays is my signal (float values from -10 to 10). What I wish to do is to normalize every datum in this array to a range of values that preceed it. (i.e. the 30001st value must have the average of the preceeding 3000 values subtracted from it and then the difference must then be divided by thisvery same average (the preceeding 3000 values). My data is collected at a rate of 100Hz thus to get a normalization of the alst 30s i must use the preceeding 3000values.

As it stand this is how I've managed to make it work:

this stores the signal into the variable photosignal

photosignal = np.array(seg.analogsignals[0], ndmin=1)

now this the part I use to get the delta F/F over a moving window of 30s

normalizedphotosignal = [(uu-(np.mean(photosignal[uu-3000:uu])))/abs(np.mean(photosignal[uu-3000:uu])) for uu in photosignal[3000:]]

The following adds 3000 values to the beginning to keep the array the same length since later on i must time lock it to another list that is the same length

holder =list(range(3000))
normalizedphotosignal = holder + normalizedphotosignal

What I have noticed is that in certain files this code gives me an error because it says that the"slice" is empty and therefore it cannot create a mean.

I think maybe there is a better way to program this that could avoid this problem altogether. Or this a correct way to approach this problem?

So i tried the solution but it is quite slow and it nevertheless still gives me the "empty slice error". I went over the moving average post and found this method: 

def running_mean(x, N):
 cumsum = np.cumsum(np.insert(x, 0, 0))
 return (cumsum[N:] - cumsum[:-N]) / N

however I'm having trouble accommodating it to my desired output. namely (x-running average)/running average

Upvotes: 2

Views: 2939

Answers (2)

Mish JPFD
Mish JPFD

Reputation: 55

Allright so I finally figured it out thanks to your help and the posts you referred me to.

The calculation for my entire data (300 000 +) takes about a second!

I used the following code:

def runningmean(x,N):
    cumsum =np.cumsum(np.insert(x,0,0))
    return (cumsum[N:] -cumsum[:-N])/N

photosignal = np.array(seg.analogsignal[0], ndmin =1) 

photosignalaverage = runningmean(photosignal, 3000) 
holder = np.zeros(2999)
photosignalaverage = np.append(holder,photosignalaverage)

detalfsignal = (photosignal-photosignalaverage)/abs(photosignalaverage)

Photosignal stores my raw signal in a numpy array. Photosignalaverage uses cumsum to calculate the running average of every datapoint in photosignal. I then add the first 2999 values as 0, to maintian the same list size as my photosignal.

I then use basic numpy calculations to get my delta F/F signal.

Thank you once more for the feedback, was truly helpful!

Upvotes: 3

phase-contrast
phase-contrast

Reputation: 116

Your approach goes in the right direction. However, you made a mistake in your list comprehension: you are using uu as your index whereas uu are the elements of your input data photosignal.

You want something like this:

normalizedphotosignal2 = np.zeros((photosignal.shape[0]-3000))
for i, uu in enumerate(photosignal[3000:]):
    normalizedphotosignal2 = (uu - (np.mean(photosignal[i-3000:i]))) / abs(np.mean(photosignal[i-3000:i]))

Keep in mind that for-loops are relatively slow in python. If performance is an issue here, you could try avoiding the for loop and use numpy methods instead (e.g. have a look at Moving average or running mean).

Hope this helps.

Upvotes: 2

Related Questions