GPU-Accelerated Polyphase Channelizer: Complete Python Implementation Guide

H_k(z) = Σ(m=0 to M-1) h(kM + m) · z^(-m)

N ≈ (As - 8)/(2.285 · Δω) + 1

@cuda.jit
def data_shuffle_kernel(input_data, output_data, Q, N):
    """Reorganize data for coalesced polyphase access"""
    idx = cuda.grid(1)
    if idx < N:
        row = idx % Q
        col = idx // Q
        output_data[row, col] = input_data[idx]

import cupy as cp
import cupyx.scipy.signal as cusignal
from numba import cuda
import numpy as np

class GPUPolyphaseChannelizer:
    def __init__(self, num_channels, filter_length, sample_rate, 
                 channel_bandwidth=None, device_id=0):
        """
        GPU-accelerated polyphase channelizer
        
        Args:
            num_channels: Number of output channels (Q)
            filter_length: Length of prototype filter (N)
            sample_rate: Input sampling rate in Hz
            channel_bandwidth: Desired channel bandwidth in Hz
            device_id: GPU device ID
        """
        self.Q = num_channels
        self.N = filter_length
        self.M = filter_length // num_channels
        self.fs = sample_rate
        self.channel_bw = channel_bandwidth or (sample_rate / num_channels)
        
        # Set GPU device
        cp.cuda.Device(device_id).use()
        
        # Design prototype filter
        self.prototype_filter = self._design_prototype_filter()
        self.polyphase_filters = self._create_polyphase_decomposition()
        
        # Pre-allocate GPU memory
        self._initialize_gpu_buffers()
        
        # Setup FFT engine
        self.fft_plan = cp.fft.get_fft_plan((self.Q,), axes=(0,))
    
    def _design_prototype_filter(self):
        """Design equiripple lowpass prototype filter"""
        # Normalized cutoff frequency
        cutoff = self.channel_bw / self.fs
        
        # Design filter using Parks-McClellan algorithm
        bands = [0, cutoff * 0.8, cutoff * 1.2, 0.5]
        desired = [1, 0]
        weights = [1, 10]  # Emphasize stopband attenuation
        
        coeffs = cusignal.remez(self.N, bands, desired, weights, Hz=1.0)
        
        # Apply window for improved stopband
        window = cusignal.get_window('kaiser', self.N, beta=8.0)
        return coeffs * window
    
    def _create_polyphase_decomposition(self):
        """Decompose prototype filter into polyphase components"""
        # Reshape filter into Q x M matrix
        polyphase = cp.zeros((self.Q, self.M), dtype=cp.float32)
        
        for q in range(self.Q):
            polyphase[q, :] = self.prototype_filter[q::self.Q]
        
        return polyphase
    
    def _initialize_gpu_buffers(self):
        """Pre-allocate GPU memory for processing"""
        # Maximum expected input size
        max_samples = 1024 * 1024
        
        self.input_buffer = cp.zeros(max_samples, dtype=cp.complex64)
        self.shuffled_buffer = cp.zeros((self.Q, max_samples // self.Q), 
                                      dtype=cp.complex64)
        self.filtered_buffer = cp.zeros((self.Q, max_samples // self.Q), 
                                      dtype=cp.complex64)
        
        # Copy filter coefficients to constant memory
        self.d_filters = cp.asarray(self.polyphase_filters)

# Polyphase filtering kernel with shared memory optimization
polyphase_filter_kernel = cp.RawKernel(r'''
extern "C" __global__
void polyphase_filter(
    const float2* __restrict__ input,
    const float* __restrict__ filters,
    float2* __restrict__ output,
    int M, int Q, int samples_per_channel)
{
    // Shared memory for filter coefficients
    extern __shared__ float shared_filters[];
    
    int channel = blockIdx.y;
    int sample = blockIdx.x * blockDim.x + threadIdx.x;
    int tid = threadIdx.x;
    
    // Cooperatively load filter coefficients
    if (tid < M) {
        shared_filters[tid] = filters[channel * M + tid];
    }
    __syncthreads();
    
    if (sample < samples_per_channel && channel < Q) {
        float2 sum = make_float2(0.0f, 0.0f);
        
        // Perform convolution using FMA operations
        #pragma unroll 4
        for (int m = 0; m < M; m++) {
            if (sample + m < samples_per_channel) {
                float2 x = input[channel * samples_per_channel + sample + m];
                float h = shared_filters[m];
                
                sum.x = fmaf(x.x, h, sum.x);
                sum.y = fmaf(x.y, h, sum.y);
            }
        }
        
        output[channel * samples_per_channel + sample] = sum;
    }
}
''', 'polyphase_filter')

class GPUPolyphaseChannelizer(GPUPolyphaseChannelizer):
    def process(self, input_signal):
        """
        Process input signal through channelizer
        
        Args:
            input_signal: Complex input samples (numpy or cupy array)
            
        Returns:
            Channelized output signals (Q channels)
        """
        # Ensure GPU array
        if isinstance(input_signal, np.ndarray):
            gpu_input = cp.asarray(input_signal, dtype=cp.complex64)
        else:
            gpu_input = input_signal
        
        n_samples = len(gpu_input)
        samples_per_channel = n_samples // self.Q
        
        # Step 1: Data shuffling for coalesced access
        self._shuffle_data(gpu_input, n_samples)
        
        # Step 2: Polyphase filtering
        threads_per_block = (256,)
        blocks_x = (samples_per_channel + 255) // 256
        blocks_per_grid = (blocks_x, self.Q)
        shared_mem_size = self.M * 4  # float32 coefficients
        
        polyphase_filter_kernel(
            blocks_per_grid, threads_per_block,
            (self.shuffled_buffer, self.d_filters, self.filtered_buffer,
             self.M, self.Q, samples_per_channel),
            shared_mem=shared_mem_size
        )
        
        # Step 3: FFT for channelization
        channelized = cp.fft.fft(self.filtered_buffer[:, :samples_per_channel], 
                                axis=0)
        
        return channelized
    
    def _shuffle_data(self, input_data, n_samples):
        """Shuffle data for polyphase decomposition"""
        samples_per_channel = n_samples // self.Q
        
        # Use CuPy's advanced indexing for efficiency
        for q in range(self.Q):
            self.shuffled_buffer[q, :samples_per_channel] = \
                input_data[q::self.Q][:samples_per_channel]

class StreamingPolyphaseChannelizer(GPUPolyphaseChannelizer):
    def __init__(self, *args, num_streams=2, **kwargs):
        super().__init__(*args, **kwargs)
        
        # Initialize CUDA streams for pipelining
        self.streams = [cp.cuda.Stream() for _ in range(num_streams)]
        self.current_stream = 0
        
        # Double buffering for continuous processing
        self.processing_buffers = [
            {
                'input': cp.zeros(8192, dtype=cp.complex64),
                'output': cp.zeros((self.Q, 8192 // self.Q), dtype=cp.complex64)
            }
            for _ in range(2)
        ]
    
    def process_streaming(self, data_generator):
        """
        Process continuous data stream
        
        Args:
            data_generator: Iterator yielding input data chunks
            
        Yields:
            Channelized output for each chunk
        """
        for chunk in data_generator:
            stream = self.streams[self.current_stream]
            buffers = self.processing_buffers[self.current_stream]
            
            with stream:
                # Asynchronous memory transfer
                buffers['input'].set(chunk, stream=stream)
                
                # Process chunk
                channelized = self.process(buffers['input'])
                
                # Copy result to output buffer
                buffers['output'][:] = channelized
            
            # Switch to next stream/buffer
            prev_stream = stream
            self.current_stream = (self.current_stream + 1) % len(self.streams)
            
            # Wait for previous processing to complete
            prev_stream.synchronize()
            
            yield buffers['output'].get(stream=stream)

# For CUDA 11.x
pip install cupy-cuda11x

# For CUDA 12.x
pip install cupy-cuda12x

# Additional dependencies
pip install numba scipy matplotlib

# Configure CuPy for optimal performance
import cupy as cp

# Enable memory pooling
mempool = cp.get_default_memory_pool()
pinned_mempool = cp.get_default_pinned_memory_pool()

# Set memory limits (e.g., 4GB)
mempool.set_limit(size=4 * 1024**3)

# Enable cuDNN (if available)
cp.cudnn.set_max_workspace_size(256 * 1024 * 1024)

# Configure for specific GPU
with cp.cuda.Device(0):
    # Get device properties
    props = cp.cuda.runtime.getDeviceProperties(0)
    print(f"GPU: {props['name'].decode()}")
    print(f"Compute Capability: {props['major']}.{props['minor']}")
    print(f"Global Memory: {props['totalGlobalMem'] / 1e9:.1f} GB")

class OptimizedMemoryManager:
    def __init__(self, max_channels=4096, max_samples=1048576):
        self.max_channels = max_channels
        self.max_samples = max_samples
        
        # Pre-allocate maximum expected memory
        self.memory_pool = {
            'input': cp.zeros(max_samples, dtype=cp.complex64),
            'shuffled': cp.zeros((max_channels, max_samples // max_channels), 
                               dtype=cp.complex64),
            'filtered': cp.zeros((max_channels, max_samples // max_channels), 
                               dtype=cp.complex64),
            'output': cp.zeros((max_channels, max_samples // max_channels), 
                             dtype=cp.complex64)
        }
        
        # Pinned memory for fast host-device transfers
        self.pinned_buffers = {
            'host_input': cp.cuda.alloc_pinned_memory(
                max_samples * np.complex64().itemsize),
            'host_output': cp.cuda.alloc_pinned_memory(
                max_channels * (max_samples // max_channels) * np.complex64().itemsize)
        }
    
    def get_buffers(self, num_channels, num_samples):
        """Get appropriately sized buffers from pool"""
        samples_per_channel = num_samples // num_channels
        
        return {
            'input': self.memory_pool['input'][:num_samples],
            'shuffled': self.memory_pool['shuffled'][:num_channels, :samples_per_channel],
            'filtered': self.memory_pool['filtered'][:num_channels, :samples_per_channel],
            'output': self.memory_pool['output'][:num_channels, :samples_per_channel]
        }

def optimize_data_transfers(channelizer, use_pinned_memory=True):
    """Configure channelizer for optimal data transfers"""
    
    if use_pinned_memory:
        # Allocate pinned memory for zero-copy access
        input_pinned = cp.cuda.alloc_pinned_memory(
            channelizer.max_input_size * cp.complex64().itemsize)
        
        # Create memory-mapped array
        channelizer.pinned_input = np.frombuffer(
            input_pinned, dtype=np.complex64).reshape(-1)
    
    # Configure unified memory for simplified management
    cp.cuda.MemPool(cp.cuda.malloc_managed).use()
    
    return channelizer

class ChannelizerBenchmark:
    def __init__(self, channelizer):
        self.channelizer = channelizer
        self.results = {}
    
    def benchmark_throughput(self, input_sizes=[1024, 8192, 65536, 524288]):
        """Measure throughput for various input sizes"""
        
        for size in input_sizes:
            # Generate test data
            test_data = cp.random.randn(size, dtype=cp.float32).astype(cp.complex64)
            
            # Warm-up
            for _ in range(5):
                self.channelizer.process(test_data)
            
            # Benchmark
            cp.cuda.Stream.null.synchronize()
            start_events = [cp.cuda.Event() for _ in range(10)]
            end_events = [cp.cuda.Event() for _ in range(10)]
            
            for i in range(10):
                start_events[i].record()
                output = self.channelizer.process(test_data)
                end_events[i].record()
            
            # Calculate statistics
            cp.cuda.Stream.null.synchronize()
            times = [cp.cuda.get_elapsed_time(start_events[i], end_events[i]) 
                    for i in range(10)]
            
            avg_time = np.mean(times) / 1000  # Convert to seconds
            throughput_msps = size / (avg_time * 1e6)
            
            self.results[size] = {
                'throughput_msps': throughput_msps,
                'latency_ms': avg_time * 1000,
                'samples_per_second': size / avg_time
            }
        
        return self.results
    
    def validate_accuracy(self, reference_implementation):
        """Validate against reference implementation"""
        test_sizes = [1024, 8192]
        
        for size in test_sizes:
            # Generate test signal with known frequency content
            freqs = [1e6, 2.5e6, 4e6]  # Hz
            test_signal = np.zeros(size, dtype=np.complex64)
            
            for f in freqs:
                t = np.arange(size) / self.channelizer.fs
                test_signal += np.exp(2j * np.pi * f * t)
            
            # Process with both implementations
            gpu_output = self.channelizer.process(cp.asarray(test_signal))
            ref_output = reference_implementation(test_signal)
            
            # Calculate error metrics
            mse = np.mean(np.abs(cp.asnumpy(gpu_output) - ref_output)**2)
            relative_error = mse / np.mean(np.abs(ref_output)**2)
            
            print(f"Size {size}: Relative Error = {relative_error:.2e}")
            assert relative_error < 1e-6, "Accuracy validation failed"

def profile_gpu_kernels(channelizer, input_data):
    """Detailed profiling of GPU kernels"""
    
    # Create NVTX ranges for profiling
    import cupy.cuda.nvtx as nvtx
    
    with nvtx.RangePush("data_shuffle"):
        channelizer._shuffle_data(input_data, len(input_data))
    
    with nvtx.RangePush("polyphase_filter"):
        # Run filter kernel
        pass  # Kernel execution here
    
    with nvtx.RangePush("fft"):
        output = cp.fft.fft(channelizer.filtered_buffer, axis=0)
    
    # Run with Nsight Systems:
    # nsys profile --stats=true python your_script.py

import pytest

class TestPolyphaseChannelizer:
    @pytest.fixture
    def channelizer(self):
        return GPUPolyphaseChannelizer(
            num_channels=16,
            filter_length=128,
            sample_rate=20e6
        )
    
    def test_filter_response(self, channelizer):
        """Test frequency response of prototype filter"""
        # Generate swept sine
        duration = 0.001  # 1ms
        t = np.arange(int(channelizer.fs * duration)) / channelizer.fs
        f0, f1 = 0, channelizer.fs / 2
        swept_sine = np.exp(2j * np.pi * (f0 + (f1 - f0) * t / duration) * t)
        
        # Process and analyze
        output = channelizer.process(cp.asarray(swept_sine))
        
        # Verify channel isolation
        for ch in range(channelizer.Q):
            channel_power = np.mean(np.abs(output[ch])**2)
            expected_band = (ch - 0.5, ch + 0.5) * channelizer.fs / channelizer.Q
            # Add assertions here
    
    def test_real_time_constraint(self, channelizer):
        """Verify real-time processing capability"""
        frame_size = 8192
        frame_duration = frame_size / channelizer.fs
        
        test_data = cp.random.randn(frame_size).astype(cp.complex64)
        
        start = time.perf_counter()
        output = channelizer.process(test_data)
        cp.cuda.Device().synchronize()
        processing_time = time.perf_counter() - start
        
        assert processing_time < frame_duration, \
            f"Real-time constraint violated: {processing_time:.3f}s > {frame_duration:.3f}s"

class GPUDebugger:
    @staticmethod
    def check_cuda_errors():
        """Check for CUDA errors after kernel execution"""
        error = cp.cuda.runtime.getLastError()
        if error != 0:
            error_name = cp.cuda.runtime.getErrorName(error)
            error_string = cp.cuda.runtime.getErrorString(error)
            raise RuntimeError(f"CUDA Error {error}: {error_name} - {error_string}")
    
    @staticmethod
    def validate_memory_access(kernel_func):
        """Decorator to validate memory access patterns"""
        def wrapper(*args, **kwargs):
            # Enable CUDA memory checking
            cp.cuda.runtime.setDevice(0)
            
            try:
                result = kernel_func(*args, **kwargs)
                cp.cuda.Device().synchronize()
                GPUDebugger.check_cuda_errors()
                return result
                
            except cp.cuda.runtime.CUDARuntimeError as e:
                print(f"Memory access error in {kernel_func.__name__}: {e}")
                # Additional debugging info
                meminfo = cp.cuda.runtime.memGetInfo()
                print(f"Free memory: {meminfo[0] / 1e9:.2f} GB")
                print(f"Total memory: {meminfo[1] / 1e9:.2f} GB")
                raise
        
        return wrapper

# Bad: Non-coalesced access
for i in range(Q):
    output[i] = input[i * stride]  # Strided access

# Good: Coalesced access with data shuffling
# First shuffle data, then access sequentially
shuffled = shuffle_for_coalescing(input)
for i in range(Q):
    output[i] = shuffled[i]  # Sequential access

def calculate_optimal_threads(num_channels, samples_per_channel):
    """Calculate optimal thread configuration"""
    
    # Ensure multiple warps per SM
    min_threads_per_block = 128
    max_threads_per_block = 1024
    
    # Optimal for coalescing
    threads_x = min(32, samples_per_channel)
    threads_y = min(max_threads_per_block // threads_x, num_channels)
    
    return (threads_x, threads_y)

from gnuradio import gr
import pmt

class gpu_channelizer_block(gr.sync_block):
    """GNU Radio block wrapper for GPU channelizer"""
    
    def __init__(self, num_channels=16, filter_length=128):
        gr.sync_block.__init__(
            self,
            name="gpu_polyphase_channelizer",
            in_sig=[np.complex64],
            out_sig=[np.complex64] * num_channels
        )
        
        self.channelizer = GPUPolyphaseChannelizer(
            num_channels=num_channels,
            filter_length=filter_length,
            sample_rate=self.sample_rate
        )
        
        self.set_output_multiple(num_channels)
    
    def work(self, input_items, output_items):
        in0 = input_items[0]
        
        # Process on GPU
        channelized = self.channelizer.process(in0)
        
        # Distribute to output streams
        for i in range(len(output_items)):
            output_items[i][:] = cp.asnumpy(channelized[i, :len(output_items[i])])
        
        return len(output_items[0])

class RealtimeChannelizerPipeline:
    def __init__(self, source_type='uhd', **kwargs):
        self.channelizer = StreamingPolyphaseChannelizer(**kwargs)
        self.source = self._create_source(source_type)
        self.running = False
    
    def _create_source(self, source_type):
        if source_type == 'uhd':
            import uhd
            usrp = uhd.usrp.MultiUSRP()
            usrp.set_rx_rate(self.channelizer.fs)
            return usrp
        elif source_type == 'rtlsdr':
            from rtlsdr import RtlSdr
            sdr = RtlSdr()
            sdr.sample_rate = self.channelizer.fs
            return sdr
    
    def start_processing(self, output_callback):
        """Start real-time processing pipeline"""
        self.running = True
        
        def data_generator():
            if hasattr(self.source, 'read_samples'):
                # RTL-SDR style
                while self.running:
                    yield self.source.read_samples(8192)
            else:
                # UHD style
                rx_stream = self.source.get_rx_stream(
                    uhd.usrp.StreamArgs('fc32', 'sc16'))
                rx_stream.issue_stream_cmd(
                    uhd.types.StreamCMD(uhd.types.StreamMode.start_cont))
                
                buffer = np.zeros(8192, dtype=np.complex64)
                while self.running:
                    rx_stream.recv(buffer, 1.0)
                    yield buffer
        
        # Process streaming data
        for channelized_output in self.channelizer.process_streaming(data_generator()):
            output_callback(channelized_output)
            
            if not self.running:
                break