ericteunis/m-thesis-introduction
1
0
Fork 0
mirror of https://gitlab.science.ru.nl/mthesis-edeboone/m-thesis-introduction.git synced 2024-12-22 11:33:32 +01:00
Code Issues Projects Releases Packages Wiki Activity

Passband calculates power not amplitude

Browse source
This commit is contained in:
Eric Teunis de Boone 2022-11-02 16:07:09 +01:00
parent 007bd7f963
commit ee0f223423
2 changed files with 46 additions and 36 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

66
fourier/04_signal_to_noise.py
View file

@ -1,4 +1,5 @@
#!/usr/bin/env python3 #!/usr/bin/env python3
# vim: fdm=indent ts=4
__doc__ = \ __doc__ = \
""" """
@ -45,31 +46,31 @@ def noisy_sine_realisation_snr(
# determine signal to noise # determine signal to noise
noise_level = bandlevel(noise, f_sample, noise_band) noise_power = bandpower(noise, f_sample, noise_band)
if cut_signal_band_from_noise_band: if cut_signal_band_from_noise_band:
lower_noise_band = passband(noise_band[0], signal_band[0]) lower_noise_band = passband(noise_band[0], signal_band[0])
upper_noise_band = passband(signal_band[1], noise_band[1]) upper_noise_band = passband(signal_band[1], noise_band[1])
noise_level = bandlevel(noise, f_sample, lower_noise_band) noise_power = bandpower(noise, f_sample, lower_noise_band)
noise_level += bandlevel(noise, f_sample, upper_noise_band) noise_power += bandpower(noise, f_sample, upper_noise_band)
signal_level = bandlevel(samples, f_sample, signal_band) signal_power = bandpower(samples, f_sample, signal_band)
snrs[j] = np.sqrt(signal_level/noise_level) snrs[j] = np.sqrt(signal_power/noise_power)
# make a nice plot showing what ranges were taken # make a nice plot showing what ranges were taken
# and the bandlevels associated with them # and the bandpowers associated with them
if return_ranges_plot and j == 0: if return_ranges_plot and j == 0:
combined_fft, freqs = ft_spectrum(samples+noise, f_sample) combined_fft, freqs = ft_spectrum(samples+noise, f_sample)
freq_scaler=1
# plot the original signal # plot the original signal
if False: if False:
_, ax = plt.subplots() _, ax = plt.subplots()
ax = plot_signal(samples+noise, sample_rate=f_sample/1e6, time_unit='us', ax=ax) ax = plot_signal(samples+noise, sample_rate=f_sample/freq_scaler, time_unit='us', ax=ax)
# plot the spectrum # plot the spectrum
if True: if True:
freq_scaler=1e6
_, axs = plot_combined_spectrum(combined_fft, freqs, freq_scaler=freq_scaler, freq_unit='MHz') _, axs = plot_combined_spectrum(combined_fft, freqs, freq_scaler=freq_scaler, freq_unit='MHz')
# indicate band ranges and frequency # indicate band ranges and frequency
@ -81,17 +82,22 @@ def noisy_sine_realisation_snr(
# indicate initial phase # indicate initial phase
axs[1].axhline(init_params[2], color='r', alpha=0.4) axs[1].axhline(init_params[2], color='r', alpha=0.4)
# plot the band levels # plot the band powers
levelax = axs[0].twinx() if False:
levelax.set_ylabel("Bandlevel") powerax = axs[0].twinx()
levelax.hlines(signal_level, noise_band[0]/freq_scaler, signal_band[1]/freq_scaler, colors=['orange']) powerax.set_ylabel("Bandpower")
levelax.hlines(noise_level, noise_band[0]/freq_scaler, noise_band[1]/freq_scaler, colors=['purple']) else:
levelax.set_ylim(bottom=0) powerax = axs[0]
powerax.hlines(np.sqrt(signal_power), noise_band[0]/freq_scaler, noise_band[1]/freq_scaler, colors=['orange'], zorder=5)
powerax.hlines(np.sqrt(noise_power), noise_band[0]/freq_scaler, noise_band[1]/freq_scaler, colors=['purple'], zorder=5)
powerax.set_ylim(bottom=0)
axs[0].legend() axs[0].legend()
# plot signal_band pass signal # plot signal_band pass signal
if False: if True:
freqs = np.fft.fftfreq(len(samples), 1/f_sample) freqs = np.fft.fftfreq(len(samples), 1/f_sample)
bandmask = bandpass_mask(freqs, band=signal_band) bandmask = bandpass_mask(freqs, band=signal_band)
fft = np.fft.fft(samples) fft = np.fft.fft(samples)
@ -99,7 +105,7 @@ def noisy_sine_realisation_snr(
bandpassed_samples = np.fft.ifft(fft) bandpassed_samples = np.fft.ifft(fft)
_, ax3 = plt.subplots() _, ax3 = plt.subplots()
ax3 = plot_signal(bandpassed_samples, sample_rate=f_sample/1e6, time_unit='us', ax=ax3) ax3 = plot_signal(bandpassed_samples, sample_rate=f_sample/freq_scaler, time_unit='us', ax=ax3)
ax3.set_title("Bandpassed Signal") ax3.set_title("Bandpassed Signal")
@ -122,10 +128,10 @@ if __name__ == "__main__":
args.fname = None args.fname = None
### ###
t_lengths = np.linspace(1e3, 5e4, 5)* 1e-9 # s t_lengths = np.linspace(1, 50, 50) # us
N = 10e1 N = 50e0
fs_sine = [33.3e6, 50e6, 73.3e6] # Hz fs_sine = [33.3, 50, 73.3] # MHz
fs_sample = [250e6, 500e6] # Hz fs_sample = [250, 500] # MHz
if False: if False:
# show t_length and fs_sample really don't care # show t_length and fs_sample really don't care
fs_iter = [ (fs_sample[0], f_sine, t_lengths) for f_sine in fs_sine ] fs_iter = [ (fs_sample[0], f_sine, t_lengths) for f_sine in fs_sine ]
@ -146,7 +152,7 @@ if __name__ == "__main__":
for i in range(int(N)): for i in range(int(N)):
delta_f = 1/t_length delta_f = 1/t_length
signal_band = passband(f_sine- 3*delta_f, f_sine + 3*delta_f) signal_band = passband(f_sine- 3*delta_f, f_sine + 3*delta_f)
noise_band = passband(30e6, 80e6) noise_band = passband(30, 80) # MHz
snrs[i], axs = noisy_sine_realisation_snr( snrs[i], axs = noisy_sine_realisation_snr(
N=1, N=1,
@ -161,7 +167,7 @@ if __name__ == "__main__":
noise_band = noise_band, noise_band = noise_band,
signal_band = signal_band, signal_band = signal_band,
return_ranges_plot=True, return_ranges_plot=False,
rng=rng, rng=rng,
) )
@ -176,14 +182,17 @@ if __name__ == "__main__":
plt.show(block=False) plt.show(block=False)
else: else:
#original code #original code
sine_amp = 1
noise_sigma = 4
my_snrs = np.zeros( (len(fs_iter), len(t_lengths), int(N)) ) my_snrs = np.zeros( (len(fs_iter), len(t_lengths), int(N)) )
for i, (f_sample, f_sine, t_lengths) in enumerate( fs_iter ): for i, (f_sample, f_sine, t_lengths) in enumerate( fs_iter ):
for k, t_length in enumerate(t_lengths): for k, t_length in enumerate(t_lengths):
return_ranges_plot = ((k==0) and True) or ( (k==(len(t_lengths)-1)) and True) and i < 1 return_ranges_plot = ((k==0) and not True) or ( (k==(len(t_lengths)-1)) and True) and i < 1
delta_f = 1/t_length delta_f = 1/t_length
signal_band = passband(f_sine- 3*delta_f, f_sine + 3*delta_f) signal_band = passband( *(f_sine + 2*delta_f*np.array([-1,1])) )
noise_band=passband(30e6, 80e6) noise_band=passband(30, 80) # MHz
my_snrs[i,k], axs = noisy_sine_realisation_snr( my_snrs[i,k], axs = noisy_sine_realisation_snr(
N=N, N=N,
@ -192,8 +201,8 @@ if __name__ == "__main__":
# signal properties # signal properties
f_sine = f_sine, f_sine = f_sine,
sine_amp = 1, sine_amp = sine_amp,
noise_sigma = 1, noise_sigma = noise_sigma,
noise_band = noise_band, noise_band = noise_band,
signal_band = signal_band, signal_band = signal_band,
@ -207,12 +216,13 @@ if __name__ == "__main__":
# plot the snrs # plot the snrs
fig, axs2 = plt.subplots() fig, axs2 = plt.subplots()
fig.basefname="signal_to_noise"
axs2.set_xlabel("$N = T*f_s$") axs2.set_xlabel("$N = T*f_s$")
axs2.set_ylabel("SNR") axs2.set_ylabel("SNR")
for i, (f_sample, f_sine, t_lengths) in enumerate(fs_iter): for i, (f_sample, f_sine, t_lengths) in enumerate(fs_iter):
# plot the means # plot the means
l = axs2.plot(t_lengths*f_sample, np.mean(my_snrs[i], axis=-1), marker='*', ls='none', label='f:{}MHz, fs:{}MHz'.format(f_sine/1e6, f_sample/1e6), markeredgecolor='black') l = axs2.plot(t_lengths*f_sample, np.mean(my_snrs[i], axis=-1), marker='*', ls='none', label='f:{}MHz, fs:{}MHz'.format(f_sine, f_sample), markeredgecolor='black')
color = l[0].get_color() color = l[0].get_color()

16
fourier/mylib/passband.py
View file

@ -16,9 +16,9 @@ class passband(namedtuple("passband", ['low', 'high'], defaults=[0, np.inf])):
def freq_mask(frequencies): def freq_mask(frequencies):
return bandpass_mask(frequencies, self) return bandpass_mask(frequencies, self)
def signal_level(samples, samplerate, normalise_bandsize=True, **ft_kwargs): def signal_power(samples, samplerate, normalise_bandsize=True, **ft_kwargs):
return bandlevel(samples, samplerate, self, normalise_bandsize, **ft_kwargs) return bandpower(samples, samplerate, self, normalise_bandsize, **ft_kwargs)
def filter_samples(samples, samplerate, **ft_kwargs): def filter_samples(samples, samplerate, **ft_kwargs):
""" """
@ -52,7 +52,7 @@ def bandpass_mask(freqs, band=passband()):
def bandsize(band = passband()): def bandsize(band = passband()):
return band[1] - band[0] return band[1] - band[0]
def bandlevel(samples, samplerate=1, band=passband(), normalise_bandsize=True, **ft_kwargs): def bandpower(samples, samplerate=1, band=passband(), normalise_bandsize=True, **ft_kwargs):
fft, freqs = ft_spectrum(samples, samplerate, **ft_kwargs) fft, freqs = ft_spectrum(samples, samplerate, **ft_kwargs)
bandmask = bandpass_mask(freqs, band=band) bandmask = bandpass_mask(freqs, band=band)
@ -62,9 +62,9 @@ def bandlevel(samples, samplerate=1, band=passband(), normalise_bandsize=True, *
else: else:
bins = 1 bins = 1
level = np.sum(np.abs(fft[bandmask])**2) power = np.sum(np.abs(fft[bandmask])**2)
return level/bins return power/bins
def signal_to_noise( samplerate, samples, noise, signal_band, noise_band=None): def signal_to_noise( samplerate, samples, noise, signal_band, noise_band=None):
if noise_band is None: if noise_band is None:
@ -73,9 +73,9 @@ def signal_to_noise( samplerate, samples, noise, signal_band, noise_band=None):
if noise is None: if noise is None:
noise = samples noise = samples
noise_level = bandlevel(noise, samplerate, noise_band) noise_power = bandpower(noise, samplerate, noise_band)
signal_level = bandlevel(samples, samplerate, signal_band) signal_power = bandpower(samples, samplerate, signal_band)
return (signal_level/noise_level)**0.5 return (signal_power/noise_power)**0.5