FFT Operations

def fft(a, n=None, axis=-1, norm=None):
    """Compute 1-D discrete Fourier Transform.
    
    Parameters:
    - a: array-like, input array
    - n: int, length of transformed axis (zero-padded or truncated)
    - axis: int, axis over which to compute FFT
    - norm: str, normalization mode ('backward', 'ortho', 'forward')
    
    Returns:
    cupy.ndarray: complex-valued FFT result
    """

def ifft(a, n=None, axis=-1, norm=None):
    """Compute 1-D inverse discrete Fourier Transform.
    
    Parameters:
    - a: array-like, input array
    - n: int, length of transformed axis
    - axis: int, axis over which to compute inverse FFT
    - norm: str, normalization mode
    
    Returns:
    cupy.ndarray: complex-valued inverse FFT result
    """

def fft2(a, s=None, axes=(-2, -1), norm=None):
    """Compute 2-D discrete Fourier Transform.
    
    Parameters:
    - a: array-like, input array
    - s: sequence of ints, shape of result along transformed axes
    - axes: sequence of ints, axes over which to compute FFT
    - norm: str, normalization mode
    
    Returns:
    cupy.ndarray: complex-valued 2D FFT result
    """

def ifft2(a, s=None, axes=(-2, -1), norm=None):
    """Compute 2-D inverse discrete Fourier Transform."""

def fftn(a, s=None, axes=None, norm=None):
    """Compute N-D discrete Fourier Transform.
    
    Parameters:
    - a: array-like, input array
    - s: sequence of ints, shape of result
    - axes: sequence of ints, axes over which to compute FFT
    - norm: str, normalization mode
    
    Returns:
    cupy.ndarray: N-dimensional FFT result
    """

def ifftn(a, s=None, axes=None, norm=None):
    """Compute N-D inverse discrete Fourier Transform."""

def rfft(a, n=None, axis=-1, norm=None):
    """Compute 1-D FFT of real input.
    
    For real input, the negative frequency terms are redundant,
    so only positive frequencies are returned.
    
    Parameters:
    - a: array-like, real-valued input array
    - n: int, length of transformed axis
    - axis: int, axis over which to compute FFT
    - norm: str, normalization mode
    
    Returns:
    cupy.ndarray: complex result with shape[axis] = n//2 + 1
    """

def irfft(a, n=None, axis=-1, norm=None):
    """Compute inverse of rfft.
    
    Parameters:
    - a: array-like, complex input from rfft
    - n: int, length of real output axis
    - axis: int, axis over which to compute inverse FFT
    - norm: str, normalization mode
    
    Returns:
    cupy.ndarray: real-valued result
    """

def rfft2(a, s=None, axes=(-2, -1), norm=None):
    """Compute 2-D FFT of real input."""

def irfft2(a, s=None, axes=(-2, -1), norm=None):
    """Compute 2-D inverse FFT to real output."""

def rfftn(a, s=None, axes=None, norm=None):
    """Compute N-D FFT of real input."""

def irfftn(a, s=None, axes=None, norm=None):
    """Compute N-D inverse FFT to real output."""

def hfft(a, n=None, axis=-1, norm=None):
    """Compute FFT of signal with Hermitian symmetry.
    
    Takes complex input with Hermitian symmetry and returns real output.
    
    Parameters:
    - a: array-like, Hermitian input array
    - n: int, length of transformed axis
    - axis: int, axis over which to compute FFT
    - norm: str, normalization mode
    
    Returns:
    cupy.ndarray: real-valued result
    """

def ihfft(a, n=None, axis=-1, norm=None):
    """Compute inverse of hfft.
    
    Takes real input and returns complex output with Hermitian symmetry.
    """

def fftfreq(n, d=1.0):
    """Return discrete Fourier Transform sample frequencies.
    
    Parameters:
    - n: int, window length
    - d: scalar, sample spacing (time between samples)
    
    Returns:
    cupy.ndarray: frequency values for each bin
    """

def rfftfreq(n, d=1.0):
    """Return sample frequencies for rfft.
    
    Parameters:
    - n: int, window length
    - d: scalar, sample spacing
    
    Returns:
    cupy.ndarray: positive frequency values
    """

def fftshift(x, axes=None):
    """Shift zero-frequency component to center of spectrum.
    
    Parameters:
    - x: array-like, input array
    - axes: int or shape tuple, axes over which to shift
    
    Returns:
    cupy.ndarray: shifted array
    """

def ifftshift(x, axes=None):
    """Inverse of fftshift.
    
    Parameters:
    - x: array-like, input array
    - axes: int or shape tuple, axes over which to shift
    
    Returns:
    cupy.ndarray: shifted array
    """

# From cupy.fft.config module
class config:
    """FFT configuration utilities for plan caching and optimization."""
    
    @staticmethod
    def get_plan_cache():
        """Get current FFT plan cache."""
    
    @staticmethod
    def set_cufft_gpus(gpus):
        """Set GPUs to use for multi-GPU FFT."""
    
    @staticmethod  
    def set_cufft_callbacks(callbacks):
        """Set cuFFT callbacks for custom processing."""

import cupy as cp
import numpy as np

# Create sample signal
n = 1024
t = cp.linspace(0, 1, n)
signal = cp.sin(2 * cp.pi * 50 * t) + 0.5 * cp.sin(2 * cp.pi * 120 * t)

# Compute FFT
fft_result = cp.fft.fft(signal)
frequencies = cp.fft.fftfreq(n, 1/n)

# Get magnitude spectrum
magnitude = cp.abs(fft_result)
print(f"Peak frequencies at bins: {cp.argmax(magnitude[:n//2])}")

import cupy as cp

# Create 2D image (or load real image data)
image = cp.random.random((512, 512))

# Compute 2D FFT
fft_image = cp.fft.fft2(image)

# Shift zero frequency to center for visualization
fft_shifted = cp.fft.fftshift(fft_image)

# Apply frequency domain filter (low-pass filter)
h, w = fft_shifted.shape
center_h, center_w = h // 2, w // 2
mask = cp.zeros((h, w))
mask[center_h-50:center_h+50, center_w-50:center_w+50] = 1

# Apply filter and inverse transform
filtered_fft = fft_shifted * mask
filtered_fft_shifted = cp.fft.ifftshift(filtered_fft)
filtered_image = cp.fft.ifft2(filtered_fft_shifted)

# Take real part (remove small imaginary components from numerical errors)
filtered_image = cp.real(filtered_image)

import cupy as cp

# Real-valued signal
real_signal = cp.random.random(10000)

# Use rfft for real input (more efficient)
rfft_result = cp.fft.rfft(real_signal)
print(f"Original length: {len(real_signal)}")
print(f"RFFT result length: {len(rfft_result)}")  # n//2 + 1

# Compute power spectral density
psd = cp.abs(rfft_result) ** 2

# Inverse transform back to real signal
reconstructed = cp.fft.irfft(rfft_result, n=len(real_signal))
print(f"Reconstruction error: {cp.mean(cp.abs(real_signal - reconstructed))}")

import cupy as cp

# Process multiple signals simultaneously
batch_size = 100
signal_length = 1024
signals = cp.random.random((batch_size, signal_length))

# FFT along last axis (each row is a signal)
batch_fft = cp.fft.fft(signals, axis=1)

# Process along first axis (each column is a signal)
column_fft = cp.fft.fft(signals, axis=0)

# 2D FFT treating each signal as 2D
signals_2d = signals.reshape(10, 10, signal_length)
fft_2d = cp.fft.fft2(signals_2d, axes=(0, 1))

import cupy as cp

# Convolution using FFT (faster for large kernels)
def fft_convolve(signal, kernel):
    """Convolution using FFT for efficiency."""
    n = len(signal) + len(kernel) - 1
    
    # Pad both to same length
    signal_padded = cp.pad(signal, (0, n - len(signal)))
    kernel_padded = cp.pad(kernel, (0, n - len(kernel)))
    
    # FFT, multiply, inverse FFT
    signal_fft = cp.fft.fft(signal_padded)
    kernel_fft = cp.fft.fft(kernel_padded)
    result_fft = signal_fft * kernel_fft
    result = cp.fft.ifft(result_fft)
    
    return cp.real(result)

# Example usage
signal = cp.random.random(10000)
gaussian_kernel = cp.exp(-cp.linspace(-3, 3, 101)**2)
convolved = fft_convolve(signal, gaussian_kernel)

# Cross-correlation using FFT
def fft_xcorr(x, y):
    """Cross-correlation using FFT."""
    n = len(x) + len(y) - 1
    x_padded = cp.pad(x, (0, n - len(x)))
    y_padded = cp.pad(y, (0, n - len(y)))
    
    X = cp.fft.fft(x_padded)
    Y = cp.fft.fft(y_padded)
    xcorr = cp.fft.ifft(X * cp.conj(Y))
    
    return cp.real(xcorr)

# Find time delay between two signals
signal1 = cp.random.random(1000)
signal2 = cp.roll(signal1, 50) + 0.1 * cp.random.random(1000)  # Delayed + noise
xcorr = fft_xcorr(signal1, signal2)
delay = cp.argmax(cp.abs(xcorr)) - len(signal1) + 1
print(f"Detected delay: {delay} samples")

import cupy as cp

def power_spectral_density(signal, fs=1.0, nperseg=256):
    """Compute power spectral density using Welch's method."""
    # Simple implementation - overlap segments and average
    n_overlap = nperseg // 2
    segments = []
    
    for i in range(0, len(signal) - nperseg, nperseg - n_overlap):
        segment = signal[i:i+nperseg]
        # Apply Hanning window
        window = 0.5 * (1 - cp.cos(2 * cp.pi * cp.arange(nperseg) / (nperseg - 1)))
        windowed = segment * window
        segments.append(windowed)
    
    if segments:
        segments = cp.stack(segments)
        # FFT each segment
        fft_segments = cp.fft.rfft(segments, axis=1)
        # Compute power and average
        power = cp.mean(cp.abs(fft_segments) ** 2, axis=0)
        freqs = cp.fft.rfftfreq(nperseg, 1/fs)
        return freqs, power
    else:
        return None, None

# Example usage
fs = 1000  # Sampling frequency
t = cp.arange(0, 10, 1/fs)
signal = cp.sin(2*cp.pi*50*t) + 0.5*cp.sin(2*cp.pi*120*t) + 0.2*cp.random.random(len(t))

freqs, psd = power_spectral_density(signal, fs)
if psd is not None:
    peak_freq_idx = cp.argmax(psd)
    peak_freq = freqs[peak_freq_idx]
    print(f"Peak frequency: {peak_freq:.1f} Hz")