FFT/Spectrogram analysis of ABF data

[swharden]

Use numpy/matplotlib to generate FFTs/spectrograms from sample ABF data

Related: https://github.com/swharden/FftSharp

[swharden]

@saglag you might find this helpful. Note that I added getAllXs() and getAllYs() just now and intend to publish an updated package later tonight. so you'll have to pip install --upgrade pyabf to get version 2.3.5

import matplotlib.pyplot as plt
import numpy as np
from typing import Tuple

def getFft(values: np.ndarray, sampleRate: int, dB: bool = False, zeroDC: bool = True) -> Tuple[np.ndarray, np.ndarray]:
    """Return FFT power and frequency for the given values

    Parameters:
        values: evenly-spaced values from the continuous signal
        sampleRate: number of signal values per second (Hz)
        dB: if True the output will be scaled to Decibel units
        zeroDC: if True the lowest frequency (0 Hz corresponding to DC offset) will be zeroed

    Returns:
        A tuple of power spectrum densities and their frequencies

    Note:
        Decibel conversion uses a power of 10 (suitable for power, not amplitude)
        https://dspillustrations.com/pages/posts/misc/decibel-conversion-factor-10-or-factor-20.html
    """

    fft = abs(np.fft.fft(values)/len(values))

    if dB:
        fft = 20 * np.log10(fft)

    if zeroDC:
        fft[0] = 0

    fftFreq = np.fft.fftfreq(len(fft), 1.0 / sampleRate)
    fft = fft[:len(fft)//2]
    fftFreq = fftFreq[:len(fftFreq)//2]
    return (fft, fftFreq)

if __name__ == "__main__":
    abfPath = pathlib.Path(PATH_ABFS).joinpath("17o05027_ic_ramp.abf")
    abf = pyabf.ABF(abfPath)
    fft, freqs = getFft(abf.getAllYs(), abf.sampleRate)

    plt.figure(figsize=(8, 4))
    plt.plot(abf.getAllXs(), abf.getAllYs())
    plt.title(f"Signal: {abf.abfID}.abf")
    plt.ylabel("Potential (mV)")
    plt.xlabel("Time (s)")
    plt.grid(alpha=.5, ls='--')
    plt.margins(0, .1)
    plt.tight_layout()

    plt.figure(figsize=(8, 4))
    plt.plot(freqs, fft)
    plt.title(f"Power Spectrum: {abf.abfID}.abf")
    plt.ylabel("Power (dB)")
    plt.xlabel("Frequency (Hz)")
    plt.grid(alpha=.5, ls='--')
    plt.axis([0, 100, 0, None])
    plt.tight_layout()
    
    plt.show()

[swharden]

@saglag following-up here's a spectrogram example. The docs have a lot more details about windowing and stepping https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.spectrogram.html

It doesn't make sense to run a spectrogram on this ABF because it's such a short recording (just a few seconds) but I don't have a longer one. If you want to email me a long recording, I'd be interested to see what that looks like 👍

import matplotlib.pyplot as plt
from scipy import signal

if __name__ == "__main__":
    abfPath = pathlib.Path(PATH_ABFS).joinpath("17o05027_ic_ramp.abf")
    abf = pyabf.ABF(abfPath)
    f, t, Sxx = signal.spectrogram(abf.getAllYs(), abf.sampleRate)
    plt.pcolormesh(t, f, Sxx)
    plt.ylabel('Frequency [Hz]')
    plt.xlabel('Time [sec]')
    plt.axis([None, None, 0, 500])
    plt.tight_layout()
    plt.show()

[saglag]

2021_09_27_0009.abf <https://drive.google.com/file/d/1AmXtlYk5wVssZVs8MceE2oBe4iHeXiB6/view?usp=drive_web> attached =D

…

On Tue, Oct 5, 2021 at 12:37 AM Scott W Harden ***@***.***> wrote: @saglag <https://github.com/saglag> following-up here's a spectrogram example. The docs have a lot more details about windowing and stepping https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.spectrogram.html It doesn't make sense to run a spectrogram on this ABF because it's such a short recording (just a few seconds) but I don't have a longer one. If you want to email me a long recording, I'd be interested to see what that looks like 👍 import matplotlib.pyplot as plt from scipy import signal if __name__ == "__main__": abfPath = pathlib.Path(PATH_ABFS).joinpath("17o05027_ic_ramp.abf") abf = pyabf.ABF(abfPath) f, t, Sxx = signal.spectrogram(abf.getAllYs(), abf.sampleRate) plt.pcolormesh(t, f, Sxx) plt.ylabel('Frequency [Hz]') plt.xlabel('Time [sec]') plt.axis([None, None, 0, 500]) plt.tight_layout() plt.show() [image: image] <https://user-images.githubusercontent.com/4165489/135960756-44cda2ac-25e7-49d3-a05d-b99017c25be6.png> — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#124 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AHW53WO3IFLTCEYE2TSI5GDUFJ6HXANCNFSM5FFI3VSQ> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.

[swharden]

Thanks! but... I don't think email attachments come through on GitHub issue pages 😅 #124

Feel free to email me directly swharden@gmail.com or use the web interface to leave a comment and drag/drop an ABF in the comment box 👍

[swharden]

My bad - I now see it's a link to a google drive folder

[saglag]

The file is unfortunately large, 1 hour at 10kHz sampling.

[swharden]

The file is unfortunately large, 1 hour at 10kHz sampling.

If we want to do spectral analysis 0-100 Hz we can downsample it to a 200 Hz sample rate without losing any frequency data 😎

A 1 hr recording of 16-bit 32-bit data at 200 Hz should be around 1.5 MB 3 MB

In ClampFit try analyze > data reduction > and decimate by whatever it takes to get to 200 Hz
(I think a sampling interval of 5000 µs?)

EDIT: after thinking about it some more I think you should use Python to do this - bin the data and take its mean (rather than decimating just taking every Nth point)

[swharden]

I poked at this a little more this morning. Here's the file it produced:

• 200hz.zip

This will be used as input to create a high quality spectrogram

Filter and Down-Sample with SciPy

• The light blue curve is the original data
• The black curve is the result of a 5-pole low-pass Bessel filter with a 100 Hz cutoff
• The orange curve is the down-sampled (200 Hz) filtered curve ideal for FFT/spectrogram analysis

NOTE: this respects the rule of thumb to sample at least 2x the low-pass filter cutoff frequency

import pyabf
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

def lowpass(data, cutoff, samplerate, order):
    normal_cutoff = cutoff / samplerate
    b, a = signal.bessel(order, normal_cutoff, btype='low', analog=False)
    y = signal.filtfilt(b, a, data)
    return y

abf = pyabf.ABF("2021_09_27_0009.abf")
ys = abf.getAllYs()
xs = np.arange(len(ys)) / abf.sampleRate
plt.plot(xs, ys, alpha=.2, label="original")

ys2 = lowpass(ys, 100, abf.sampleRate, 5)
plt.plot(xs, ys2, 'k', label="filtered")

newSampleRate = 200
reductionFactor = abf.sampleRate // newSampleRate
ys3 = ys2[::reductionFactor]
xs3 = np.arange(len(ys3)) / newSampleRate
plt.plot(xs3, ys3, label="downsampled")

plt.axis([1850, 1860, None, None])
plt.show()

np.save("200hz.npy", ys3)

[swharden]

Refer to https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.spectrogram.html

• increase vertical resolution by increasing nfft (power of 2)
• increase vertical sensitivity by modifying nperseg (power of 2)
• increase horizontal resolution by changing noverlap

Notice I used that low-pass-filtered down-sampled (200 Hz) .npy file as my data input, but there's no reason this wouldn't work on the raw 20 kHz data (it would just be slower to calculate).

Here's a comparison of the FFT frequencies at two different time points in the experiment:

Here's a spectrogram showing spectral power over time:

Here's the code that generates these two images:

import numpy as np
from scipy import signal
from scipy.fft import fftshift
import matplotlib.pyplot as plt

ys = np.load("200hz.npy")
sampleRate = 200

f, t, Sxx = signal.spectrogram(
    ys, sampleRate, nfft=pow(2, 13), nperseg=pow(2, 9))

# compare FFTs at 2 time points
timeA = 30  # minutes
timeB = 90  # minutes
fftA = Sxx[int(timeA / t[1])]
fftB = Sxx[int(timeB / t[1])]
fftA[0] = 0  # zero the DC component
fftB[0] = 0  # zero the DC component
fftFreq = np.arange(len(fftA)) / sampleRate

# low-pass filter the FFT to make peaks more broad
b, a = signal.butter(3, 0.04)  # tweak for smoothness
fftA = signal.filtfilt(b, a, fftA)
fftB = signal.filtfilt(b, a, fftB)

plt.plot(fftFreq, fftA, alpha=.8, label=f"minute {timeA}")
plt.plot(fftFreq, fftB, alpha=.8, label=f"minute {timeB}")
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power')
plt.legend()
plt.grid(alpha=.5, ls='--')
plt.axis([None, None, 0, None])
plt.tight_layout()
plt.show()

# show FFTs over time as a spectrogram
im = plt.pcolormesh(t / 60, fftshift(f),
                    fftshift(Sxx, axes=0), shading='gouraud')
plt.colorbar(im)
plt.ylabel('Frequency (Hz)')
plt.xlabel('Time (minutes)')
plt.axis([None, None, 0, 10])
plt.tight_layout()
plt.show()

This is a little more nuanced than ideal for the pyabf tutorial so I'll probably just leave this information here (it's googlable).

These images and code snippets should be enough to help people get off the ground with ABF frequency analysis 👍

Thanks the suggestion @saglag! It was interesting to apply signal analysis to ephys data in this way 🚀