|
| 1 | +"""Spectral analysis helper module.""" |
| 2 | + |
| 3 | +import numpy as np |
| 4 | + |
| 5 | + |
| 6 | +def db_to_pow(value, places=3): |
| 7 | + """Convert dBW to W.""" |
| 8 | + if isinstance(value, np.ndarray): |
| 9 | + return np.round(10 ** (0.1 * value), places) |
| 10 | + return round(10 ** (0.1 * value), places) |
| 11 | + |
| 12 | + |
| 13 | +def dBW(value, places=1): |
| 14 | + """Convert to dBW.""" |
| 15 | + if isinstance(value, np.ndarray): |
| 16 | + return np.round(10 * np.log10(value), places) |
| 17 | + return round(10 * np.log10(value), places) |
| 18 | + |
| 19 | + |
| 20 | +def enob(sndr, places=1): |
| 21 | + """Return ENOB for given SNDR.""" |
| 22 | + return round((sndr - 1.76) / 6.02, places) |
| 23 | + |
| 24 | + |
| 25 | +def sndr_sfdr(spectrum, freq, fs, nfft, leak, full_scale=0): |
| 26 | + """Get SNDR and SFDR.""" |
| 27 | + |
| 28 | + # Zero the DC bin |
| 29 | + for i in range(0, leak + 1): |
| 30 | + spectrum[i] = 0 |
| 31 | + bin_sig = np.argmax(spectrum) |
| 32 | + psig = sum(spectrum[i] for i in range(bin_sig - leak, bin_sig + leak + 1)) |
| 33 | + spectrum_n = spectrum.copy() |
| 34 | + spectrum_n[bin_sig] = 0 |
| 35 | + fbin = fs / nfft |
| 36 | + |
| 37 | + for i in range(bin_sig - leak, bin_sig + leak + 1): |
| 38 | + spectrum_n[i] = 0 |
| 39 | + |
| 40 | + bin_spur = np.argmax(spectrum_n) |
| 41 | + pspur = spectrum[bin_spur] |
| 42 | + |
| 43 | + noise_power = sum(spectrum_n) |
| 44 | + noise_floor = 2 * noise_power / nfft |
| 45 | + |
| 46 | + stats = {} |
| 47 | + |
| 48 | + stats["sig"] = { |
| 49 | + "freq": freq[bin_sig], |
| 50 | + "bin": bin_sig, |
| 51 | + "power": psig, |
| 52 | + "dB": dBW(psig), |
| 53 | + "dBFS": round(dBW(psig) - full_scale, 1), |
| 54 | + } |
| 55 | + |
| 56 | + stats["spur"] = { |
| 57 | + "freq": freq[bin_spur], |
| 58 | + "bin": bin_spur, |
| 59 | + "power": pspur, |
| 60 | + "dB": dBW(pspur), |
| 61 | + "dBFS": round(dBW(pspur) - full_scale, 1), |
| 62 | + } |
| 63 | + stats["noise"] = { |
| 64 | + "floor": noise_floor, |
| 65 | + "power": noise_power, |
| 66 | + "rms": np.sqrt(noise_power), |
| 67 | + "dBHz": round(dBW(noise_floor, 3) - full_scale, 1), |
| 68 | + "NSD": round(dBW(noise_floor, 3) - full_scale - 2 * dBW(fbin, 3), 1), |
| 69 | + } |
| 70 | + stats["sndr"] = { |
| 71 | + "dBc": dBW(psig / noise_power), |
| 72 | + "dBFS": round(full_scale - dBW(noise_power), 1), |
| 73 | + } |
| 74 | + stats["sfdr"] = { |
| 75 | + "dBc": dBW(psig / pspur), |
| 76 | + "dBFS": round(full_scale - dBW(pspur), 1), |
| 77 | + } |
| 78 | + stats["enob"] = {"bits": enob(stats["sndr"]["dBFS"])} |
| 79 | + |
| 80 | + return stats |
| 81 | + |
| 82 | + |
| 83 | +def find_harmonics(spectrum, freq, nfft, bin_sig, psig, harms=5, leak=20, fscale=1e6): |
| 84 | + """Get the harmonic contents of the data.""" |
| 85 | + harm_stats = {"harm": {}} |
| 86 | + harm_index = 2 |
| 87 | + for harm in bin_sig * np.arange(2, harms + 1): |
| 88 | + harm_stats["harm"][harm_index] = {} |
| 89 | + zone = np.floor(harm / (nfft / 2)) + 1 |
| 90 | + if zone % 2 == 0: |
| 91 | + bin_harm = int(nfft / 2 - (harm - (zone - 1) * nfft / 2)) |
| 92 | + else: |
| 93 | + bin_harm = int(harm - (zone - 1) * nfft / 2) |
| 94 | + |
| 95 | + # Make sure we pick the max bin where power is maximized; due to spectral leakage |
| 96 | + # if bin_harm == nfft/2, set to bin of 0 |
| 97 | + if bin_harm == nfft / 2: |
| 98 | + bin_harm = 0 |
| 99 | + pwr_max = spectrum[bin_harm] |
| 100 | + bin_harm_max = bin_harm |
| 101 | + for i in range(bin_harm - leak, bin_harm + leak + 1): |
| 102 | + try: |
| 103 | + pwr = spectrum[i] |
| 104 | + if pwr > pwr_max: |
| 105 | + bin_harm_max = i |
| 106 | + pwr_max = pwr |
| 107 | + except IndexError: |
| 108 | + # bin + leakage out of bounds, so stop looking |
| 109 | + break |
| 110 | + |
| 111 | + harm_stats["harm"][harm_index]["bin"] = bin_harm_max |
| 112 | + harm_stats["harm"][harm_index]["power"] = pwr_max |
| 113 | + harm_stats["harm"][harm_index]["freq"] = round(freq[bin_harm] / fscale, 1) |
| 114 | + harm_stats["harm"][harm_index]["dBc"] = dBW(pwr_max / psig) |
| 115 | + harm_stats["harm"][harm_index]["dB"] = dBW(pwr_max) |
| 116 | + |
| 117 | + harm_index = harm_index + 1 |
| 118 | + |
| 119 | + return harm_stats |
| 120 | + |
| 121 | + |
| 122 | +def get_plot_string(stats, full_scale, fs, nfft, window, xscale=1e6, fscale="MHz"): |
| 123 | + """Generate plot string from stats dict.""" |
| 124 | + |
| 125 | + plt_str = "==== FFT ====\n" |
| 126 | + plt_str += f"NFFT = {nfft}\n" |
| 127 | + plt_str += f"fbin = {round(fs/nfft / 1e3, 2)} kHz\n" |
| 128 | + plt_str += f"window = {window}\n" |
| 129 | + plt_str += "\n" |
| 130 | + plt_str += "==== Signal ====\n" |
| 131 | + plt_str += f"FullScale = {full_scale} dB\n" |
| 132 | + plt_str += f"Psig = {stats['sig']['dBFS']} dBFS ({stats['sig']['dB']} dB)\n" |
| 133 | + plt_str += f"fsig = {round(stats['sig']['freq']/xscale, 2)} {fscale}\n" |
| 134 | + plt_str += f"fsamp = {round(fs/xscale, 2)} {fscale}\n" |
| 135 | + plt_str += "\n" |
| 136 | + plt_str += "==== SNDR/SFDR ====\n" |
| 137 | + plt_str += f"ENOB = {stats['enob']['bits']} bits\n" |
| 138 | + plt_str += f"SNDR = {stats['sndr']['dBFS']} dBFS ({stats['sndr']['dBc']} dBc)\n" |
| 139 | + plt_str += f"SFDR = {stats['sfdr']['dBFS']} dBFS ({stats['sfdr']['dBc']} dBc)\n" |
| 140 | + plt_str += f"Pspur = {stats['spur']['dBFS']} dBFS\n" |
| 141 | + plt_str += f"fspur = {round(stats['spur']['freq']/xscale, 2)} {fscale}\n" |
| 142 | + plt_str += f"Noise Floor = {stats['noise']['dBHz']} dBFS\n" |
| 143 | + plt_str += f"NSD = {stats['noise']['NSD']} dBFS\n" |
| 144 | + plt_str += "\n" |
| 145 | + plt_str += "==== Harmonics ====\n" |
| 146 | + |
| 147 | + for hindex, hdata in stats["harm"].items(): |
| 148 | + plt_str += f"HD{hindex} = {round(hdata['dB'] - full_scale, 1)} dBFS @ {hdata['freq']} {fscale}\n" |
| 149 | + |
| 150 | + plt_str += "\n" |
| 151 | + |
| 152 | + return plt_str |
0 commit comments