# -*- coding: utf-8 -*- # ######### COPYRIGHT ######### # Credits # ####### # # Copyright(c) 2015-2025 # ---------------------- # # * `LabEx Archimède `_ # * `Laboratoire d'Informatique Fondamentale `_ # (now `Laboratoire d'Informatique et Systèmes `_) # * `Institut de Mathématiques de Marseille `_ # * `Université d'Aix-Marseille `_ # # This software is a port from LTFAT 2.1.0 : # Copyright (C) 2005-2025 Peter L. Soendergaard . # # Contributors # ------------ # # * Denis Arrivault # * Florent Jaillet # # Description # ----------- # # ltfatpy is a partial Python port of the # `Large Time/Frequency Analysis Toolbox `_, # a MATLAB®/Octave toolbox for working with time-frequency analysis and # synthesis. # # Version # ------- # # * ltfatpy version = 1.1.2 # * LTFAT version = 2.1.0 # # Licence # ------- # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program. If not, see . # # ######### COPYRIGHT ######### """This module contains samples of a periodized Hermite function Ported from ltfat_2.1.0/fourier/pherm.m .. moduleauthor:: Denis Arrivault """ import numpy as np from ltfatpy.comp.comp_hermite import comp_hermite from ltfatpy.comp.comp_hermite_all import comp_hermite_all from ltfatpy.sigproc.normalize import normalize def pherm(L, order, tfr=1, phase='accurate', orthtype='noorth'): """PHERM Periodized Hermite function - Usage: | ``g, = pherm(L,order)`` | ``g, = pherm(L,order,tfr)`` | ``g, D = pherm(...)`` - Input parameters: :param int L: Length of vector. :param order: Order of Hermite function. :type order: scalar or numpy.ndarray :param float tfr: ratio between time and frequency support. 1 by default :param str phase: 'accurate' or 'fast' (see below) :param str orthtype: 'noorth', 'polar' or 'qr' (see below). - Output parameters: :returns: ``(g, D)`` :rtype: tuple :var numpy.ndarray g: The periodized Hermite function :var numpy.ndarray D: The eigenvalues of the Discrete Fourier Transform corresponding to the Hermite functions. ``pherm(L,order,tfr)`` computes samples of a periodized Hermite function of order **order**. **order** is counted from 0, so the zero'th order Hermite function is the Gaussian. The parameter **tfr** determines the ratio between the effective support of **g** and the effective support of the DFT of **g**. If :math:`tfr>1` then **g** has a wider support than the DFT of **g**. ``pherm(L,order)`` does the same setting :math:`tfr=1`. If **order** is a vector, ``pherm`` will return a matrix, where each column is a Hermite function with the corresponding **order**. ``g, D = pherm(...)`` also returns the eigenvalues **D** of the Discrete Fourier Transform corresponding to the Hermite functions. The returned functions are eigenvectors of the DFT. The Hermite functions are orthogonal to all other Hermite functions with a different eigenvalue, but eigenvectors with the same eigenvalue are not orthogonal (but see the flags below). **phase** can take the following values: 'accurate' By default it uses a numerically very accurate that computes each Hermite function individually. This is the default 'fast' Use a less accurate algorithm that calculates all the Hermite up to a given order at once. **orthtype** can take the following values: 'noorth' orthonormalization of the Hermite functions. This is the default. 'polar' Orthonormalization of the Hermite functions using the polar decomposition orthonormalization method. 'qr' Orthonormalization of the Hermite functions using the Gram-Schmidt orthonormalization method (usign ``qr``). If you just need to compute a single Hermite function, there is no speed difference between the **accurate** and **fast** algorithm. - Examples: The following plot shows the spectrograms of 4 Hermite functions of length 200 with order 1, 10, 100, and 190::: >>> import numpy as np >>> import matplotlib.pyplot as plt >>> from ltfatpy import sgram >>> plt.close('all') >>> _ = plt.figure() >>> _ = plt.subplot(221) >>> _ = sgram(pherm(200, 1)[0], nf=True, tc=True, normalization='lin', ... colorbar=False) >>> _ = plt.subplot(2,2,2) >>> _ = sgram(pherm(200, 10)[0], nf=True, tc=True, normalization='lin', ... colorbar=False) >>> _ = plt.subplot(2,2,3) >>> _ = sgram(pherm(200, 100)[0], nf=True, tc=True, ... normalization='lin', colorbar=False) >>> _ = plt.subplot(2,2,4) >>> _ = sgram(pherm(200, 190)[0], nf=True, tc=True, ... normalization='lin', colorbar=False) >>> plt.show() .. image:: images/pherm.png :width: 700px :alt: spectrograms :align: center .. seealso:: :func:`~ltfatpy.fourier.pgauss.pgauss`, :func:`~ltfatpy.fourier.psech.psech` """ if not np.isscalar(L) or isinstance(L, str): raise TypeError("L must be a scalar") if not isinstance(L, int): raise TypeError('L must be an integer') # Parse tfr and order. if (not np.isscalar(tfr)): raise TypeError('tfr must be a scalar or vector') if np.isscalar(order) and not isinstance(order, str): W = 1 order = np.array([order]) elif isinstance(order, np.ndarray): order = order.reshape(-1).copy() W = order.shape[0] else: raise TypeError('order must be a scalar or vector') # Calculate W windows. if 'accurate' in phase: # Calculate W windows. g = np.zeros((L, W)) for w in range(W): thisorder = order[w] safe = get_safe(thisorder) # Outside the interval [-safe,safe] # then H(thisorder) is numerically zero. nk = int(np.ceil(safe/np.sqrt(L/np.sqrt(tfr)))) sqrtl = np.sqrt(L) lr = np.arange(L) for k in range(-nk, nk+1): xval = (lr/sqrtl - k*sqrtl) / np.sqrt(tfr) g[:, w] = g[:, w] + comp_hermite(thisorder, np.sqrt(2*np.pi)*xval) else: highestorder = np.max(order) safe = get_safe(highestorder) # Outside the interval [-safe,safe] # then H(thisorder) is numerically zero. nk = int(np.ceil(safe/np.sqrt(L/np.sqrt(tfr)))) g = np.zeros((L, highestorder+1)) sqrtl = np.sqrt(L) lr = np.arange(L) for k in range(-nk, nk+1): xval = (lr/sqrtl - k*sqrtl)/np.sqrt(tfr) g = g + comp_hermite_all(highestorder+1, np.sqrt(2*np.pi)*xval) g = g[:, order] if 'polar' in orthtype: # Orthonormalize within each of the 4 eigenspaces for ii in range(4): subidx = ((order % 4) == ii) gsub = g[:, subidx] if gsub.size: U, _, V = np.linalg.svd(gsub, full_matrices=False) gsub = np.dot(U, V) else: gsub = np.asarray([]) g[:, subidx] = gsub if 'qr' in orthtype: # Orthonormalize within each of the 4 eigenspaces for ii in range(4): subidx = ((order % 4) == ii) gsub = g[:, subidx] if gsub.size: Q, _ = np.linalg.qr(gsub, mode='reduced') else: Q = np.asarray([]) g[:, subidx] = Q if 'noorth' in orthtype: # Just normalize it, no orthonormalization g = normalize(g)[0] # set up the eigenvalues D = np.exp(-1j*order*np.pi/2) if W == 1: g = g.squeeze() return(g, D) def get_safe(order): # These numbers have been computed numerically. if order <= 6: safe = 4 else: if order <= 18: safe = 5 else: if order <= 31: safe = 6 else: if order <= 46: safe = 7 else: # Anything else, use a high number. safe = 12 return safe if __name__ == '__main__': # pragma: no cover import doctest doctest.testmod()