psychopy / psychopy

For running psychology and neuroscience experiments
http://www.psychopy.org
GNU General Public License v3.0
1.64k stars 897 forks source link

[Feature Request]: Clarification of Units for cutoff in Butterworth Filters #5875

Closed altunenes closed 9 months ago

altunenes commented 10 months ago

Description

Hi and ty for psychopy :)

I've been diving into the psychopy/visual/filters.py file and I've found the Butterworth filters to be incredibly useful. However, I might have missed something, but it seems that the units for the cutoff parameter aren't explicitly documented in this section. While the cutoffvalue appears to be normalized between 0 and 1, I wasn't able to find a clear indication or functionality that translates this into more standardized units like cycles per pixel or cycles per degree of visual angle.

If I've overlooked any existing documentation or functionality regarding this, I'd greatly appreciate being pointed in the right direction.

Purpose

No response

Is your feature request related to a problem?

A little bit of detailed documentation for the cutoff value (or maybe a function for unit conversion) could be beneficial for researchers and reviewers when they are working with SFs.

Describe the solution you would like

If my understanding is correct, the cutoff value might be multiplied by 0.5 to convert it into cycles per pixel. I base this assumption on the fact that the maximum frequency in a 128x128 image, for example, is 64 cycles, or 0.5 cycles per pixel. Thus, a cutoff of 0.2 might correspond to 0.2Γ—0.5=0.1 cycles per pixel (?)

So If my assumptions hold true (?), a function within filters.py for unit conversion would be a valuable addition. This could even extend to converting cycles per pixel to cycles per degree, requiring parameters like screen/img dimensions and viewing distance etc etc

Describe alternatives that you have considered

No response

Additional context

I can be a volunteer if help is needed, thankssss :=)

peircej commented 10 months ago

That code was added many years ago by a contributor but, from reading it, I don't believe it knows anything about real-world units coming from PsychoPy so its unit is presumably either pixels or fractions of the filter. I think the issue is that, when that function is used it isn't clear where its output will be applied. It's just the mathematical function

altunenes commented 10 months ago

thanks @peircej.

In fact, Spatial Frequency Filtering example seems a comprehensive guide on this topic I mean it shows us how can we use these filters for applying SFs ❀️ . The example begins by taking an image and converting it into the frequency domain using the fft. Once the image is in the frequency domain, its amplitude spectrum is calculated using np.fft.fftshift(np.abs(img_freq)) (where img_freq is the fft transformed image).

The code uses the psychopy.filters.butter2d_lp function and specifies both the size of the image (raw_img.shape) and the cutoff frequency (0.05 in this case). Additionally, the 'order' of the filter is defined, which controls the steepness of the filter curve (order of Butterworth)

lp_filt = psychopy.filters.butter2d_lp(
    size=raw_img.shape,
    cutoff=0.05,
    n=10
)

Once the Butterworth filter is configured, it's then applied to the frequency-space representation of the image using element-wise multiplication: img_filt = np.fft.fftshift(img_freq) * lp_filt. Following this, the amplitude of the filtered image is taken, and certain transformations are applied to visualize the result/image. :=)

While the example is incredibly useful, one crucial aspect seems to be missing: it doesn't clearly state the units of the cutoff parameter in the documentation. This is particularly important for those who want to understand how the cutoffs translate into cycles per pixel (so we can transform it into c/degree). πŸ₯²

thanks πŸ™‚

peircej commented 10 months ago

To reiterate, I don't know the answer but I'm pretty sure it is either 1/pixels, or 1/filter_width. Maybe you could deduce which of those is correct from your own empirical tests and add to the documentation.

altunenes commented 10 months ago

Thanks πŸ™‚ In my first spare time, I'll run some tests and loop back here once I've got it figured out. πŸ§‘β€πŸŽ“

altunenes commented 9 months ago

TL/DR: After diving into several signal processing textbooks and resources, I can confirm that cylce/pixel=cutoff*0.5 in psychopy.

The Butterworth filter can be formulated in various ways, but This formulation, although slightly different in terms of variable naming, is essentially identical.

For instance, consider this formula from a widely referenced textbook [1]:

image

B low is the filter response (Butterworth) D(u,v) represents the distance from the origin in the frequency domain. D_0 is the cutoff frequency n is the order of the filter.

Another formulate I would like to pick from a vision journal [2] )

image

In both formulas, despite their slight differences, their roles and the overall expression remain analogous.

when we examine the psychopy's formula:

f = 1 / (1.0 + (radius/cutoff)**(2 * n))

which can we transform it into latex:

image

note: R(u,v) corresponds to the distance from the center of the 2D frequency plane,

but what does the "cutoff" frequency mean in these formulas?

In the formulas and practical applications of Butterworth filters, the cutoff frequency is expressed as a ratio, effectively nullifying the physical units, be it samples per second (Hertz) or samples per meter, leaving us with cycles per sample. When we talk about image processing, these 'samples' correspond to pixels. Here, the Nyquist Sampling Theorem plays a pivotal role, stipulating that a value of 0.5 aligns with half a cycle per sample, thereby defining the upper Nyquist limit [5].

In simpler terms, a cutoff of 1.0 implies that the Butterworth filter will permit all frequency components up to the Nyquist limit to pass through without attenuation, thereby rendering an unfiltered signal. Consequently, the maximum Nyquist frequency consistently stands at 0.5 cycles/pixel.

So, in the context of our Butterworth filter in psychopy, the cutoff parameter, normalized between 0 and 1, denotes the relative cutoff frequency of the filter. A cutoff of 1.0 indicates an unimpeded passage for all frequency components up to the Nyquist limit, essentially reproducing an unaltered signal. This conceptualization brings forth the relation:

cycle/pixel = cutoff Γ— 0.5

In bellow , I made a 1D signal and used a 1D version of the Butterworth filter from Psychopy

import numpy as np
import matplotlib.pyplot as plt
from numpy.fft import fft, ifft

# simple sinusoidal signal..
x = np.linspace(0, 4 * np.pi, 500) 
y = np.sin(x)                       

# Butterworth filter in 1D.. (in psychopy we have 2D)
def butterworth_1d(cutoff, n=3):
    radius = np.linspace(-0.5, 0.5, len(x))
    f = 1 / (1.0 + (radius/cutoff)**(2*n))
    return f

cutoffs = [1.0, 0.5, 0.1]
filtered_signals = []

for cutoff in cutoffs:
    filter_response = butterworth_1d(cutoff)
    filtered_y = np.real(ifft(fft(y) * filter_response))
    filtered_signals.append(filtered_y)

plt.figure(figsize=(12, 6))
plt.plot(x, y, label='Original Signal', color='black', linewidth=2)
colors = ['red', 'green', 'blue']
for cutoff, filtered_y, color in zip(cutoffs, filtered_signals, colors):
    plt.plot(x, filtered_y, label=f'Cutoff = {cutoff}', color=color)
plt.title('Effect of Butterworth Filter on a Sinusoidal Signal')
plt.xlabel('x')
plt.ylabel('Amplitude')
plt.legend()
plt.grid(True)
plt.show()

image

when you look at the 1.0 cutoff value which is the maximum cutoff value it's the same as the original signal. This means it allows all of the original signal frequencies to pass through, so it looks exactly like the original signal. lower cutoff signal gets more and more smoothed as the cutoff gets smaller because fewer frequencies are allowed to pass through the filter...

I hope this description will be helpful for those who want to use Butterworth filters in psychopy...

[1] - Sonka, M., Hlavac, V.,, Boyle, R. (1998). Image Processing: Analysis and Machine Vision. CL-Engineering. ISBN: 053495393X [2] - Kwon, M., & Legge, G. E. (2011). Spatial-frequency cutoff requirements for pattern recognition in central and peripheral vision. Vision research, 51(18), 1995–2007. https://doi.org/10.1016/j.visres.2011.06.020 [3] Chapter 10, The Essential Physics of Medical Imaging,Bushberg Radiation Detection and Measurement [4] Carlsson, G. A. and Ljungberg, M. (2021). Basic atomic and nuclear physics. Handbook of Nuclear Medicine and Molecular Imaging for Physicists, 15-37. https://doi.org/10.1201/9780429489556-2 [5] Oppenheim, A. V., Schafer, R. W.,, Buck, J. R. (1999). Discrete-Time Signal Processing. Prentice-hall Englewood Cliffs.

some websites (about Nyquist sampling theorem) or <- https://www.kth.se/social/files/542d1224f27654321376a4b4/Compendium.Imaging.Physics.pdf or <- http://www.people.vcu.edu/~mhcrosthwait/clrs322/SPECTimagefiltering.htm