Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Remove phase for FFT normalizations 'power' and 'psd' #557

Merged
merged 3 commits into from
Apr 12, 2024
Merged

Remove phase for FFT normalizations 'power' and 'psd' #557

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
bc91efd
Update fft.py
Mar 12, 2024
43faa92
Update test_fft.py
Mar 12, 2024
6928b97
Update test_fft.py
Mar 12, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 2 additions & 8 deletions pyfar/dsp/fft.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -229,23 +229,17 @@ def normalization(spec, n_samples, sampling_rate, fft_norm='none',
else:
# Equation 12 in Ahrens et al. 2020
norm /= np.sum(window)**2
# the phase is kept for being able to switch between normalizations
# altoug the power spectrum does usually not have phase information,
# i.e., spec = np.abs(spec)**2
if not inverse:
spec *= np.abs(spec)
spec = np.abs(spec)**2
elif fft_norm == 'psd':
if window is None:
# Equation 6 in Ahrens et al. 2020
norm /= (n_samples * sampling_rate)
else:
# Equation 13 in Ahrens et al. 2020
norm /= (np.sum(np.asarray(window)**2) * sampling_rate)
# the phase is kept for being able to switch between normalizations
# altoug the power spectrum does usually not have phase information,
# i.e., spec = np.abs(spec)**2
if not inverse:
spec *= np.abs(spec)
spec = np.abs(spec)**2
elif fft_norm != 'unitary':
raise ValueError(("norm type must be 'unitary', 'amplitude', 'rms', "
f"'power', or 'psd' but is '{fft_norm}'"))
Expand Down
32 changes: 24 additions & 8 deletions tests/test_fft.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -139,7 +139,7 @@ def test_normalization_none(impulse_stub):
def test_normalization_single_sided_single_channel_even_samples():
# single sided test spectrum
v = 1/3 + 1/3j
vsq = v * np.abs(v)
vsq = np.abs(v)**2
spec_single = np.array([v, v, v])
# valid number of samples of time signal corresponding to spec_single
N = 4 # time signal with even number of samples
Expand Down Expand Up @@ -173,13 +173,17 @@ def test_normalization_single_sided_single_channel_even_samples():
print(f"Assesing normalization: '{normalization}' (inverse)")
spec_out_inv = fft.normalization(spec_out, N, fs,
normalization, inverse=True)
npt.assert_allclose(spec_out_inv, spec_single, atol=1e-15)
if normalization in ['power', 'psd']:
# Inverse only resembles magnitude for power norms
npt.assert_allclose(spec_out_inv, np.abs(spec_single), atol=1e-15)
else:
npt.assert_allclose(spec_out_inv, spec_single, atol=1e-15)


def test_normalization_single_sided_single_channel_odd_samples():
# single sided test spectrum
v = 1/3 + 1/3j
vsq = v * np.abs(v)
vsq = np.abs(v)**2
spec_single = np.array([v, v, v])
# valid number of samples of time signal corresponding to spec_single
N = 5 # time signal with even number of samples
Expand Down Expand Up @@ -213,13 +217,17 @@ def test_normalization_single_sided_single_channel_odd_samples():
print(f"Assesing normalization: '{normalization}' (inverse)")
spec_out_inv = fft.normalization(spec_out, N, fs,
normalization, inverse=True)
npt.assert_allclose(spec_out_inv, spec_single, atol=1e-15)
if normalization in ['power', 'psd']:
# Inverse only resembles magnitude for power norms
npt.assert_allclose(spec_out_inv, np.abs(spec_single), atol=1e-15)
else:
npt.assert_allclose(spec_out_inv, spec_single, atol=1e-15)


def test_normalization_both_sided_single_channel():
# single sided test spectrum
v = 1/3 + 1/3j
vsq = v * np.abs(v)
vsq = np.abs(v)**2
spec_single = np.array([v, v, v])
# valid number of samples of time signal corresponding to spec_single
N = 3 # time signal with even number of samples
Expand Down Expand Up @@ -251,13 +259,17 @@ def test_normalization_both_sided_single_channel():
spec_out_inv = fft.normalization(spec_out, N, fs,
normalization, inverse=True,
single_sided=False)
npt.assert_allclose(spec_out_inv, spec_single, atol=1e-15)
if normalization in ['power', 'psd']:
# Inverse only resembles magnitude for power norms
npt.assert_allclose(spec_out_inv, np.abs(spec_single), atol=1e-15)
else:
npt.assert_allclose(spec_out_inv, spec_single, atol=1e-15)


def test_normalization_single_sided_multi_channel_even_samples():
# single sided test spectrum
v = 1/3 + 1/3j
vsq = v * np.abs(v)
vsq = np.abs(v)**2
tile = (4, 2, 1)
spec_single = np.tile(np.array([v, v, v]), tile)
# valid number of samples of time signal corresponding to spec_single
Expand Down Expand Up @@ -293,7 +305,11 @@ def test_normalization_single_sided_multi_channel_even_samples():
print(f"Assesing normalization: '{normalization}' (inverse)")
spec_out_inv = fft.normalization(spec_out, N, fs,
normalization, inverse=True)
npt.assert_allclose(spec_out_inv, spec_single, atol=1e-15)
if normalization in ['power', 'psd']:
# Inverse only resembles magnitude for power norms
npt.assert_allclose(spec_out_inv, np.abs(spec_single), atol=1e-15)
else:
npt.assert_allclose(spec_out_inv, spec_single, atol=1e-15)


def test_normalization_with_window():
Expand Down