|
| 1 | +#!/usr/bin/env python3 |
| 2 | +# -*- coding: utf-8 -*- |
| 3 | + |
| 4 | +""" |
| 5 | +Some basic spectral analysis. |
| 6 | +
|
| 7 | +References: |
| 8 | + - Analog Devices MT-003 TUTORIAL |
| 9 | + "Understand SINAD, ENOB, SNR, THD, THD + N, and SFDR so You Don't Get Lost in the Noise Floor" |
| 10 | + - National Instruments Application Note 041 |
| 11 | + "The Fundamentals of FFT-Based Signal Analysis and Measurement" |
| 12 | +""" |
| 13 | + |
| 14 | +import numpy as np |
| 15 | +import matplotlib |
| 16 | +import matplotlib.pyplot as plt |
| 17 | + |
| 18 | + |
| 19 | +def db2amp(db): |
| 20 | + """Decibels to amplitutde ratio""" |
| 21 | + return 10 ** (0.05 * db) |
| 22 | + |
| 23 | + |
| 24 | +def amp2db(a): |
| 25 | + """Amplitutde ratio to decibels""" |
| 26 | + return 20 * np.log10(a) |
| 27 | + |
| 28 | + |
| 29 | +def db2pow(db): |
| 30 | + """Decibels to power ratio""" |
| 31 | + return 10 ** (0.1 * db) |
| 32 | + |
| 33 | + |
| 34 | +def pow2db(p): |
| 35 | + """Power ratio to decibels""" |
| 36 | + return 10 * np.log10(p) |
| 37 | + |
| 38 | + |
| 39 | +def enob(sinad): |
| 40 | + """Calculate ENOB from SINAD""" |
| 41 | + return (sinad - 1.76) / 6.02 |
| 42 | + |
| 43 | + |
| 44 | +def snr_theor(n): |
| 45 | + """Theoretical SNR of an ideal n-bit ADC in dB""" |
| 46 | + return 6.02 * n + 1.76 |
| 47 | + |
| 48 | + |
| 49 | +def noise_floor(snr, m): |
| 50 | + """Noise floor of the m-point FFT in dB""" |
| 51 | + return -snr - 10 * np.log10(m / 2) |
| 52 | + |
| 53 | + |
| 54 | +def harmonics(psp, fft_n, ref_pow, sample_freq, leak=20, n=5, window='hanning'): |
| 55 | + """Obtain first n harmonics properties from power spectrum""" |
| 56 | + # Coherence Gain and Noise Power Bandwidth for different windows |
| 57 | + win_params = {'uniform': {'cg': 1.0, 'npb': 1.0}, |
| 58 | + 'hanning': {'cg': 0.5, 'npb': 1.5}, |
| 59 | + 'hamming': {'cg': 0.54, 'npb': 1.36}, |
| 60 | + 'blackman': {'cg': 0.42, 'npb': 1.73}}[window] |
| 61 | + fft_n = len(psp) * 2 # one side spectrum provided |
| 62 | + df = sample_freq / fft_n |
| 63 | + # calculate fundamental frequency |
| 64 | + fund_bin = np.argmax(psp) |
| 65 | + fund_freq = (np.sum([psp[i] * i * df for i in range(fund_bin - leak, fund_bin + leak + 1)]) / |
| 66 | + np.sum(psp[fund_bin - leak: fund_bin + leak + 1])) |
| 67 | + # calculate harmonics info |
| 68 | + h = [] |
| 69 | + for i in range(1, n + 1): |
| 70 | + h_i = {'num': i} |
| 71 | + zone_freq = (fund_freq * i) % sample_freq |
| 72 | + h_i['freq'] = sample_freq - zone_freq if zone_freq >= (sample_freq / 2) else zone_freq |
| 73 | + h_i['central_bin'] = int(h_i['freq'] / df) |
| 74 | + h_i['bins'] = np.array(range(h_i['central_bin'] - leak, h_i['central_bin'] + leak + 1)) |
| 75 | + h_i['pow'] = ((1 / win_params['cg']) ** 2) * np.sum(psp[h_i['bins']]) / win_params['npb'] |
| 76 | + h_i['vrms'] = np.sqrt(h_i['pow']) |
| 77 | + if i == 1: |
| 78 | + h_i['db'] = '%.2f dBFS' % pow2db(h_i['pow'] / ref_pow) |
| 79 | + else: |
| 80 | + h_i['db'] = '%.2f dBc' % pow2db(h_i['pow'] / h[0]['pow']) |
| 81 | + h += [h_i] |
| 82 | + return h |
| 83 | + |
| 84 | + |
| 85 | +def signal_noise(psp, harmonics): |
| 86 | + """Obtain different signal+noise characteristics from spectrum""" |
| 87 | + # noise + distortion power |
| 88 | + nd_psp = np.copy(psp) |
| 89 | + nd_psp[harmonics[0]['bins']] = 0 # remove main harmonic |
| 90 | + nd_psp[0] = 0 # remove dc |
| 91 | + nd_pow = sum(nd_psp) |
| 92 | + # noise power |
| 93 | + n_psp = np.copy(psp) |
| 94 | + for h in harmonics: |
| 95 | + n_psp[h['bins']] = 0 # remove all harmonics |
| 96 | + n_psp[0] = 0 # remove dc |
| 97 | + n_pow = sum(n_psp) |
| 98 | + # distortion power |
| 99 | + d_pow = np.sum([h['pow'] for h in harmonics]) - harmonics[0]['pow'] |
| 100 | + # calculate results |
| 101 | + sinad = pow2db(harmonics[0]['pow'] / nd_pow) |
| 102 | + thd = pow2db(harmonics[0]['pow'] / d_pow) |
| 103 | + snr = pow2db(harmonics[0]['pow'] / n_pow) |
| 104 | + sfdr = pow2db(max(nd_psp) / harmonics[0]['pow']) |
| 105 | + return sinad, thd, snr, sfdr |
| 106 | + |
| 107 | + |
| 108 | +def analyze(sig, adc_bits, adc_vref, adc_freq, window='hanning'): |
| 109 | + """Do spectral analysis for ADC samples""" |
| 110 | + # Calculate some useful parameters |
| 111 | + sig_vpeak_max = adc_vref / 2 |
| 112 | + sig_vrms_max = sig_vpeak_max / np.sqrt(2) |
| 113 | + sig_pow_max = sig_vrms_max ** 2 |
| 114 | + ref_pow = sig_pow_max |
| 115 | + adc_prd = 1 / adc_freq |
| 116 | + adc_quants = 2 ** adc_bits |
| 117 | + dv = adc_vref / adc_quants |
| 118 | + sig_n = len(sig) |
| 119 | + dt = 1 / adc_freq |
| 120 | + fft_n = sig_n |
| 121 | + df = adc_freq / fft_n |
| 122 | + win_coef = {'uniform': np.ones(sig_n), |
| 123 | + 'hanning': np.hanning(sig_n)}[window] |
| 124 | + sp_leak = 20 # spectru leak bins |
| 125 | + h_n = 5 # harmonics number |
| 126 | + |
| 127 | + # Convert samples to voltage |
| 128 | + sig_v = sig * dv |
| 129 | + |
| 130 | + # Remove DC and apply window |
| 131 | + sig_dc = np.mean(sig_v) |
| 132 | + sig_windowed = (sig_v - sig_dc) * win_coef |
| 133 | + |
| 134 | + # Calculate one-side amplitude spectrum (Vrms) |
| 135 | + asp = np.sqrt(2) * np.abs(np.fft.rfft(sig_windowed)) / sig_n |
| 136 | + |
| 137 | + # Calculate one-side power spectrum (Vrms^2) |
| 138 | + psp = np.power(asp, 2) |
| 139 | + psp_db = pow2db(psp / ref_pow) |
| 140 | + |
| 141 | + # Calculate harmonics |
| 142 | + h = harmonics(psp=psp, fft_n=fft_n, ref_pow=ref_pow, sample_freq=adc_freq, leak=sp_leak, n=h_n, window=window) |
| 143 | + |
| 144 | + # Input signal parameters (based on 1st harmonic) |
| 145 | + sig_pow = h[0]['pow'] |
| 146 | + sig_vrms = h[0]['vrms'] |
| 147 | + sig_vpeak = sig_vrms * np.sqrt(2) |
| 148 | + sig_freq = h[0]['freq'] |
| 149 | + sig_prd = 1 / sig_freq |
| 150 | + |
| 151 | + # Calculate SINAD, THD, SNR, SFDR |
| 152 | + adc_sinad, adc_thd, adc_snr, adc_sfdr = signal_noise(psp, h) |
| 153 | + |
| 154 | + # Calculate ENOB |
| 155 | + # sinad correction to normalize ENOB to full-scale regardless of input signal amplitude |
| 156 | + adc_enob = enob(adc_sinad + pow2db(ref_pow / sig_pow)) |
| 157 | + |
| 158 | + # Calculate Noise Floor |
| 159 | + adc_noise_floor = noise_floor(adc_snr, fft_n) |
| 160 | + |
| 161 | + # Create plots |
| 162 | + fig = plt.figure(figsize=(14, 7)) |
| 163 | + gs = matplotlib.gridspec.GridSpec(2, 2, width_ratios=[3, 1]) |
| 164 | + |
| 165 | + # Time plot |
| 166 | + ax_time = plt.subplot(gs[0, 0]) |
| 167 | + ax_time_xlim = min(sig_n, int(5 * sig_prd / dt)) |
| 168 | + ax_time.plot(np.arange(0, ax_time_xlim), sig[:ax_time_xlim], color='C0') |
| 169 | + ax_time.set(ylabel='ADC code', ylim=[0, adc_quants]) |
| 170 | + ax_time.set(yticks=list(range(0, adc_quants, adc_quants // 8)) + [adc_quants - 1]) |
| 171 | + ax_time.set(xlabel='Sample', xlim=[0, ax_time_xlim - 1]) |
| 172 | + ax_time.set(xticks=range(0, ax_time_xlim, max(1, ax_time_xlim // 20))) |
| 173 | + ax_time.grid(True) |
| 174 | + ax_time_xsec = ax_time.twiny() |
| 175 | + ax_time_xsec.set(xticks=ax_time.get_xticks()) |
| 176 | + ax_time_xsec.set(xbound=ax_time.get_xbound()) |
| 177 | + ax_time_xsec.set_xticklabels(['%.02f' % (x * dt * 1e3) for x in ax_time.get_xticks()]) |
| 178 | + ax_time_xsec.set_xlabel('Time, ms') |
| 179 | + ax_time_ysec = ax_time.twinx() |
| 180 | + ax_time_ysec.set(yticks=ax_time.get_yticks()) |
| 181 | + ax_time_ysec.set(ybound=ax_time.get_ybound()) |
| 182 | + ax_time_ysec.set_yticklabels(['%.02f' % (x * dv) for x in ax_time.get_yticks()]) |
| 183 | + ax_time_ysec.set_ylabel('Voltage, V') |
| 184 | + |
| 185 | + # Frequency plot |
| 186 | + ax_freq = plt.subplot(gs[1, 0]) |
| 187 | + ax_freq.plot(np.arange(0, len(psp_db)), psp_db, color='C0', zorder=0, label="Spectrum") |
| 188 | + for h_i in h: |
| 189 | + ax_freq.text(h_i['central_bin'] + 2, psp_db[h_i['central_bin']], str(h_i['num']), |
| 190 | + va='bottom', ha='left', weight='bold') |
| 191 | + ax_freq.plot(h_i['bins'], psp_db[h_i['bins']], color='C4') |
| 192 | + ax_freq.plot(0, 0, color='C4', label="Harmonics") |
| 193 | + ax_freq.set(ylabel='dB', ylim=[-150, 10]) |
| 194 | + ax_freq.set(xlabel='Sample', xlim=[0, fft_n / 2]) |
| 195 | + ax_freq.set(xticks=list(range(0, fft_n // 2, fft_n // 32)) + [fft_n // 2 - 1]) |
| 196 | + ax_freq.grid(True) |
| 197 | + ax_freq.legend(loc="lower right", ncol=3) |
| 198 | + ax_freq_sec = ax_freq.twiny() |
| 199 | + ax_freq_sec.set_xticks(ax_freq.get_xticks()) |
| 200 | + ax_freq_sec.set_xbound(ax_freq.get_xbound()) |
| 201 | + ax_freq_sec.set_xticklabels(['%.02f' % (x * df * 1e-3) for x in ax_freq.get_xticks()]) |
| 202 | + ax_freq_sec.set_xlabel('Frequency, kHz') |
| 203 | + |
| 204 | + # Information plot |
| 205 | + ax_info = plt.subplot(gs[:, 1]) |
| 206 | + ax_info.set(xlim=[0, 10], xticks=[], ylim=[0, 10], yticks=[]) |
| 207 | + harmonics_str = '\n'.join(['%d%s @ %-10s : %s' % (h_i['num'], ['st', 'nd', 'rd', 'th', 'th'][h_i['num'] - 1], |
| 208 | + '%0.3f kHz' % (h_i['freq'] * 1e-3), |
| 209 | + h_i['db']) for h_i in h]) |
| 210 | + ax_info_str = """ |
| 211 | +========= FFT ========== |
| 212 | +Points : {fft_n} |
| 213 | +Freq. resolution : {fft_res:.4} Hz |
| 214 | +Window : {fft_window} |
| 215 | +
|
| 216 | +======= Harmonics ====== |
| 217 | +{harmonics_str} |
| 218 | +
|
| 219 | +===== Input signal ===== |
| 220 | +Frequency : {sig_freq:.4} kHz |
| 221 | +Amplitude (Vpeak): {sig_vpeak:.4} V |
| 222 | +DC offset : {sig_dc:.4} V |
| 223 | +
|
| 224 | +========= ADC ========== |
| 225 | +Sampling freq. : {adc_freq:.4} kHz |
| 226 | +Sampling period : {adc_prd:.4} us |
| 227 | +Reference volt. : {adc_vref:.4} V |
| 228 | +Bits : {adc_bits} bits |
| 229 | +Quants : {adc_quants} |
| 230 | +Quant : {adc_quant:.4} mV |
| 231 | +SNR : {adc_snr:.4} dB |
| 232 | +SINAD : {adc_sinad:.4} dB |
| 233 | +THD : {adc_thd:.4} dB |
| 234 | +ENOB : {adc_enob:.4} bits |
| 235 | +SFDR : {adc_sfdr:.4} dBc |
| 236 | +Noise floor : {adc_nfloor:.4} dBFS |
| 237 | +""".format(fft_n=fft_n, |
| 238 | + fft_res=df, |
| 239 | + fft_window=window, |
| 240 | + harmonics_str=harmonics_str, |
| 241 | + sig_freq=sig_freq * 1e-3, |
| 242 | + sig_vpeak=sig_vpeak, |
| 243 | + sig_dc=sig_dc, |
| 244 | + adc_freq=adc_freq * 1e-3, |
| 245 | + adc_prd=adc_prd * 1e6, |
| 246 | + adc_vref=adc_vref, |
| 247 | + adc_bits=adc_bits, |
| 248 | + adc_quants=adc_quants, |
| 249 | + adc_quant=dv * 1e3, |
| 250 | + adc_snr=adc_snr, |
| 251 | + adc_thd=adc_thd, |
| 252 | + adc_sinad=adc_sinad, |
| 253 | + adc_enob=adc_enob, |
| 254 | + adc_sfdr=adc_sfdr, |
| 255 | + adc_nfloor=adc_noise_floor) |
| 256 | + ax_info.text(1, 9.5, ax_info_str, va='top', ha='left', family='monospace') |
| 257 | + |
| 258 | + # General plotting settings |
| 259 | + plt.tight_layout() |
| 260 | + plt.style.use('bmh') |
| 261 | + |
| 262 | + # Show the result |
| 263 | + plt.show() |
0 commit comments