def getFFT(signal,T):
N = len(signal)
xf = np.linspace(0.0, 1.0/(2.0*float(T)), N/2)
w = blackman(N)
ffts = []
sfft = fft(np.array(signal)*w)
return xf,sfft
python类blackman()的实例源码
def getFFT(signal,T):
N = len(signal)
xf = np.linspace(0.0, 1.0/(2.0*float(T)), N/2)
w = blackman(N)
ffts = []
sfft = fft(np.array(signal)*w)
return xf,sfft
def plot_spectrum(sys):
"""
Plot the High Harmonic Generation spectrum
"""
# Power spectrum emitted is calculated using the Larmor formula
# (https://en.wikipedia.org/wiki/Larmor_formula)
# which says that the power emitted is proportional to the square of the acceleration
# i.e., the RHS of the second Ehrenfest theorem
N = len(sys.P_average_RHS)
k = np.arange(N)
# frequency range
omegas = (k - N / 2) * np.pi / (0.5 * sys.t)
# spectra of the
spectrum = np.abs(
# used windows fourier transform to calculate the spectra
# rhttp://docs.scipy.org/doc/scipy/reference/tutorial/fftpack.html
fftpack.fft((-1) ** k * blackman(N) * sys.P_average_RHS)
) ** 2
spectrum /= spectrum.max()
plt.semilogy(omegas / sys.omega_laser, spectrum)
plt.ylabel('spectrum (arbitrary units)')
plt.xlabel('frequency / $\\omega_L$')
plt.xlim([0, 45.])
plt.ylim([1e-15, 1.])
def get_ft_windows():
""" Retrieve the available windows to be applied on signal data before FT.
@return: dict with keys being the window name and items being again a dict
containing the actual function and the normalization factor to
calculate correctly the amplitude spectrum in the Fourier Transform
To find out the amplitude normalization factor check either the scipy
implementation on
https://github.com/scipy/scipy/blob/v0.15.1/scipy/signal/windows.py#L336
or just perform a sum of the window (oscillating parts of the window should
be averaged out and constant offset factor will remain):
MM=1000000 # choose a big number
print(sum(signal.hanning(MM))/MM)
"""
win = {'none': {'func': np.ones, 'ampl_norm': 1.0},
'hamming': {'func': signal.hamming, 'ampl_norm': 1.0/0.54},
'hann': {'func': signal.hann, 'ampl_norm': 1.0/0.5},
'blackman': {'func': signal.blackman, 'ampl_norm': 1.0/0.42},
'triang': {'func': signal.triang, 'ampl_norm': 1.0/0.5},
'flattop': {'func': signal.flattop, 'ampl_norm': 1.0/0.2156},
'bartlett': {'func': signal.bartlett, 'ampl_norm': 1.0/0.5},
'parzen': {'func': signal.parzen, 'ampl_norm': 1.0/0.375},
'bohman': {'func': signal.bohman, 'ampl_norm': 1.0/0.4052847},
'blackmanharris': {'func': signal.blackmanharris, 'ampl_norm': 1.0/0.35875},
'nuttall': {'func': signal.nuttall, 'ampl_norm': 1.0/0.3635819},
'barthann': {'func': signal.barthann, 'ampl_norm': 1.0/0.5}}
return win
def stft(time_signal, time_dim=None, size=1024, shift=256,
window=signal.blackman, fading=True, window_length=None):
"""
Calculates the short time Fourier transformation of a multi channel multi
speaker time signal. It is able to add additional zeros for fade-in and
fade out and should yield an STFT signal which allows perfect
reconstruction.
:param time_signal: multi channel time signal.
:param time_dim: Scalar dim of time.
Default: None means the biggest dimension
:param size: Scalar FFT-size.
:param shift: Scalar FFT-shift. Typically shift is a fraction of size.
:param window: Window function handle.
:param fading: Pads the signal with zeros for better reconstruction.
:param window_length: Sometimes one desires to use a shorter window than
the fft size. In that case, the window is padded with zeros.
The default is to use the fft-size as a window size.
:return: Single channel complex STFT signal
with dimensions frames times size/2+1.
"""
if time_dim is None:
time_dim = np.argmax(time_signal.shape)
# Pad with zeros to have enough samples for the window function to fade.
if fading:
pad = [(0, 0)] * time_signal.ndim
pad[time_dim] = [size - shift, size - shift]
time_signal = np.pad(time_signal, pad, mode='constant')
# Pad with trailing zeros, to have an integral number of frames.
frames = _samples_to_stft_frames(time_signal.shape[time_dim], size, shift)
samples = _stft_frames_to_samples(frames, size, shift)
pad = [(0, 0)] * time_signal.ndim
pad[time_dim] = [0, samples - time_signal.shape[time_dim]]
time_signal = np.pad(time_signal, pad, mode='constant')
if window_length is None:
window = window(size)
else:
window = window(window_length)
window = np.pad(window, (0, size - window_length), mode='constant')
time_signal_seg = segment_axis(time_signal, size,
size - shift, axis=time_dim)
letters = string.ascii_lowercase
mapping = letters[:time_signal_seg.ndim] + ',' + letters[time_dim + 1] \
+ '->' + letters[:time_signal_seg.ndim]
return rfft(np.einsum(mapping, time_signal_seg, window),
axis=time_dim + 1)
def istft(stft_signal, size=1024, shift=256,
window=signal.blackman, fading=True, window_length=None):
"""
Calculated the inverse short time Fourier transform to exactly reconstruct
the time signal.
:param stft_signal: Single channel complex STFT signal
with dimensions frames times size/2+1.
:param size: Scalar FFT-size.
:param shift: Scalar FFT-shift. Typically shift is a fraction of size.
:param window: Window function handle.
:param fading: Removes the additional padding, if done during STFT.
:param window_length: Sometimes one desires to use a shorter window than
the fft size. In that case, the window is padded with zeros.
The default is to use the fft-size as a window size.
:return: Single channel complex STFT signal
:return: Single channel time signal.
"""
assert stft_signal.shape[1] == size // 2 + 1
if window_length is None:
window = window(size)
else:
window = window(window_length)
window = np.pad(window, (0, size - window_length), mode='constant')
window = _biorthogonal_window_loopy(window, shift)
# Why? Line created by Hai, Lukas does not know, why it exists.
window *= size
time_signal = scipy.zeros(stft_signal.shape[0] * shift + size - shift)
for j, i in enumerate(range(0, len(time_signal) - size + shift, shift)):
time_signal[i:i + size] += window * np.real(irfft(stft_signal[j]))
# Compensate fade-in and fade-out
if fading:
time_signal = time_signal[
size - shift:len(time_signal) - (size - shift)]
return time_signal
def stft(time_signal, time_dim=None, size=1024, shift=256,
window=signal.blackman, fading=True, window_length=None):
"""
Calculates the short time Fourier transformation of a multi channel multi
speaker time signal. It is able to add additional zeros for fade-in and
fade out and should yield an STFT signal which allows perfect
reconstruction.
:param time_signal: multi channel time signal.
:param time_dim: Scalar dim of time.
Default: None means the biggest dimension
:param size: Scalar FFT-size.
:param shift: Scalar FFT-shift. Typically shift is a fraction of size.
:param window: Window function handle.
:param fading: Pads the signal with zeros for better reconstruction.
:param window_length: Sometimes one desires to use a shorter window than
the fft size. In that case, the window is padded with zeros.
The default is to use the fft-size as a window size.
:return: Single channel complex STFT signal
with dimensions frames times size/2+1.
"""
if time_dim is None:
time_dim = np.argmax(time_signal.shape)
'''
# Pad with zeros to have enough samples for the window function to fade.
if fading:
pad = [(0, 0)] * time_signal.ndim
pad[time_dim] = [size - shift, size - shift]
time_signal = np.pad(time_signal, pad, mode='constant')
'''
# Pad with trailing zeros, to have an integral number of frames.
frames = _samples_to_stft_frames(time_signal.shape[time_dim], size, shift)
samples = _stft_frames_to_samples(frames, size, shift)
pad = [(0, 0)] * time_signal.ndim
pad[time_dim] = [0, samples - time_signal.shape[time_dim]]
time_signal = np.pad(time_signal, pad, mode='constant')
if window_length is None:
window = window(size)
else:
window = window(window_length)
window = np.pad(window, (0, size - window_length), mode='constant')
time_signal_seg = segment_axis(time_signal, size,
size - shift, axis=time_dim)
letters = string.ascii_lowercase
mapping = letters[:time_signal_seg.ndim] + ',' + letters[time_dim + 1] \
+ '->' + letters[:time_signal_seg.ndim]
return rfft(np.einsum(mapping, time_signal_seg, window),
axis=time_dim + 1)
def istft(stft_signal, signal_len, size=1024, shift=256,
window=signal.blackman, fading=True, window_length=None):
"""
Calculated the inverse short time Fourier transform to exactly reconstruct
the time signal.
:param stft_signal: Single channel complex STFT signal
with dimensions frames times size/2+1.
:param size: Scalar FFT-size.
:param shift: Scalar FFT-shift. Typically shift is a fraction of size.
:param window: Window function handle.
:param fading: Removes the additional padding, if done during STFT.
:param window_length: Sometimes one desires to use a shorter window than
the fft size. In that case, the window is padded with zeros.
The default is to use the fft-size as a window size.
:return: Single channel complex STFT signal
:return: Single channel time signal.
"""
assert stft_signal.shape[1] == size // 2 + 1
if window_length is None:
window = window(size)
else:
window = window(window_length)
window = np.pad(window, (0, size - window_length), mode='constant')
window = _biorthogonal_window_loopy(window, shift)
# Why? Line created by Hai, Lukas does not know, why it exists.
window *= size
time_signal = scipy.zeros(stft_signal.shape[0] * shift + size - shift)
for j, i in enumerate(range(0, len(time_signal) - size + shift, shift)):
time_signal[i:i + size] += window * np.real(irfft(stft_signal[j]))
time_signal = time_signal[0:signal_len]
# Compensate fade-in and fade-out
#if fading:
# time_signal = time_signal[
# size - shift:len(time_signal) - (size - shift)]
return time_signal
def stft(time_signal, time_dim=None, size=1024, shift=256,
window=signal.blackman, fading=True, window_length=None):
"""
Calculates the short time Fourier transformation of a multi channel multi
speaker time signal. It is able to add additional zeros for fade-in and
fade out and should yield an STFT signal which allows perfect
reconstruction.
:param time_signal: multi channel time signal.
:param time_dim: Scalar dim of time.
Default: None means the biggest dimension
:param size: Scalar FFT-size.
:param shift: Scalar FFT-shift. Typically shift is a fraction of size.
:param window: Window function handle.
:param fading: Pads the signal with zeros for better reconstruction.
:param window_length: Sometimes one desires to use a shorter window than
the fft size. In that case, the window is padded with zeros.
The default is to use the fft-size as a window size.
:return: Single channel complex STFT signal
with dimensions frames times size/2+1.
"""
if time_dim is None:
time_dim = np.argmax(time_signal.shape)
# Pad with zeros to have enough samples for the window function to fade.
if fading:
pad = [(0, 0)] * time_signal.ndim
pad[time_dim] = [size - shift, size - shift]
time_signal = np.pad(time_signal, pad, mode='constant')
# Pad with trailing zeros, to have an integral number of frames.
frames = _samples_to_stft_frames(time_signal.shape[time_dim], size, shift)
samples = _stft_frames_to_samples(frames, size, shift)
pad = [(0, 0)] * time_signal.ndim
pad[time_dim] = [0, samples - time_signal.shape[time_dim]]
time_signal = np.pad(time_signal, pad, mode='constant')
if window_length is None:
window = window(size)
else:
window = window(window_length)
window = np.pad(window, (0, size - window_length), mode='constant')
time_signal_seg = segment_axis(time_signal, size,
size - shift, axis=time_dim)
letters = string.ascii_lowercase
mapping = letters[:time_signal_seg.ndim] + ',' + letters[time_dim + 1] \
+ '->' + letters[:time_signal_seg.ndim]
return rfft(np.einsum(mapping, time_signal_seg, window),
axis=time_dim + 1)
def istft(stft_signal, size=1024, shift=256,
window=signal.blackman, fading=True, window_length=None):
"""
Calculated the inverse short time Fourier transform to exactly reconstruct
the time signal.
:param stft_signal: Single channel complex STFT signal
with dimensions frames times size/2+1.
:param size: Scalar FFT-size.
:param shift: Scalar FFT-shift. Typically shift is a fraction of size.
:param window: Window function handle.
:param fading: Removes the additional padding, if done during STFT.
:param window_length: Sometimes one desires to use a shorter window than
the fft size. In that case, the window is padded with zeros.
The default is to use the fft-size as a window size.
:return: Single channel complex STFT signal
:return: Single channel time signal.
"""
assert stft_signal.shape[1] == size // 2 + 1
if window_length is None:
window = window(size)
else:
window = window(window_length)
window = np.pad(window, (0, size - window_length), mode='constant')
window = _biorthogonal_window_loopy(window, shift)
# Why? Line created by Hai, Lukas does not know, why it exists.
window *= size
time_signal = scipy.zeros(stft_signal.shape[0] * shift + size - shift)
for j, i in enumerate(range(0, len(time_signal) - size + shift, shift)):
time_signal[i:i + size] += window * np.real(irfft(stft_signal[j]))
# Compensate fade-in and fade-out
if fading:
time_signal = time_signal[
size - shift:len(time_signal) - (size - shift)]
return time_signal
