Coverage for ana/mfn.py : 39%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# SPDX-FileCopyrightText: 2020 Jonathan Pieper <jpieper@stud.uni-frankfurt.de>
2#
3# SPDX-License-Identifier: GPL-3.0-or-later
5"""Handles multiple Magnetic Flux Noise (MFN) Measurements
6taken with the SR830 DAQ."""
8import logging
9import os
10from glob import glob
12import matplotlib.pyplot as plt
13import numpy as np
14import pandas as pd
15import scipy.stats
16import seaborn as sns
17from matplotlib.colors import LogNorm
19from .raw import RAW
20from .multiple import MultipleM
23def load_meas_info(file='data/mfn_info.csv'):
24 meas_info = pd.read_csv(file, index_col=0)
25 meas_info.index.name = 'Nr'
27 # Drop Measurements where an error occured.
28 meas_info = meas_info[meas_info['Error'] != True]
29 meas_info.drop('Error', axis=1, inplace=True)
30 return meas_info
32class MFN(MultipleM):
33 def __init__(self, files, **kwargs):
34 r"""Magnetic Flux Noise Measurements.
35 Class to handle and plot MFN Measurements
36 taken with the SR830 DAQ.
38 if a measurement number is given the files are
39 loaded from folder:
40 ``data/MFN/m%s/*.dat``
42 else: files must be a list of filenames:
44 .. code:: python
45 :linenos:
46 :caption: Example
48 arg = glob('data/MFN/mXXX/m[0-9.]*.dat')
49 mfn = MFN(arg)
51 :param files: list of files or measurement number.
52 :param kwargs: Different additional parameters.
53 :param int nr: Measurement number
54 (default: None).
55 :param bool calc_first: Calculating first spectrum
56 (default: True).
57 :param bool calc_second: Calculating second spectrum
58 (default: False).
59 :param bool downsample: downsample timesignal for first spectrum
60 (default: False).
61 :param int num_first_spectra: Number of spectra to average
62 (default: 64).
63 :param float timeconstant: Lock-In timeconstant to cut first spectrum
64 at corresponding frequency (see ``cut_fmax``).
65 :param float cut_fmax: Cut first spectrum at this frequency
66 (default: :math:`f_{max} = \frac{\tau}{2\pi}`).
67 :param default_cm: matplotlib color map for contour plots
68 (default: style['default_cm']).
69 :type default_cm: matplotlib.cm
70 """
71 super().__init__()
72 self.name = 'Magnetic Flux Noise'
73 self.default_cm = kwargs.get('default_cm', self.style.get('default_cm'))
74 self.style['default_cm'] = self.default_cm
76 self.info = {
77 'Nr': kwargs.get('nr'),
78 'first': kwargs.get('calc_first', True),
79 'second': kwargs.get('calc_second', False),
80 'num': kwargs.get('num_first_spectra', 64),
81 'downsample': kwargs.get('downsample', False),
82 'length': kwargs.get('length'),
83 }
85 if isinstance(files, int):
86 self.info['Nr'] = files
87 files = glob('data/MFN/m%s/*' % self.info['Nr'])
89 self.logger = logging.getLogger('MFNMeasurement (%s)' %
90 self.info['Nr'])
92 self.data, self.measurements = self.load_meas(files, **kwargs)
94 # Stop here if we have no data
95 if len(self.data) == 0:
96 return
98 # Cut too long measurements
99 if kwargs.get('equalize_length', True):
100 self.__equalize_length()
102 self.info['timeconstant'] = kwargs.get('timeconstant')
104 # Cut higher frequencies if timeconstant is known
105 fmax_rate = self.info.get('rate') / (2*np.pi) if self.info.get('rate') else 1e6
106 fmax_tc = 1 / (2 * np.pi * self.info['timeconstant']) if self.info['timeconstant'] else 1e6
107 fmax = np.min([fmax_rate, fmax_tc])
108 self.info['fmax'] = fmax
109 if fmax < 1e6 or kwargs.get('cut_fmax'):
110 self.cut_fmax(kwargs.get('cut_fmax', fmax))
112 def load_meas(self, files, **kwargs):
113 """
114 Loading a single magnetic flux noise DAQ measurement.
116 Files in measrument folder contain multiple measurements
117 with one measurement per field.
119 :param files: files to load
121 :return: DataFrame Data, dict Measurement Objects
122 :rtype: tuple
123 """
125 max_len = 0
126 meas = {}
127 all_data = pd.DataFrame({'Field': [], 'Time': [], 'Vx': [], 'Vy': []})
128 for f in files:
129 info_dict = self.get_info_from_name(f)
131 field = info_dict.get('field')
132 nr = info_dict.get('nr')
134 if not self.info.get('Nr'):
135 self.info['Nr'] = nr
137 data_df = pd.read_csv(f, sep='\t')
138 if data_df.empty:
139 continue
141 if len(data_df['Vx']) % 1024:
142 d = data_df['Vx'].iloc[:-(len(data_df['Vx']) % 1024)]
143 else:
144 d = data_df.Vx
146 max_len = len(d)
148 data = {
149 'data': d,
150 'info': {
151 'Nr': nr,
152 'field': field,
153 'rate': 1 / data_df.time.diff().mean(),
154 'length': max_len * data_df.time.diff().mean(),
155 }
156 }
158 if not self.info['length']:
159 self.info['length'] = max_len * data_df.time.diff().mean()
160 if not self.info.get('rate'):
161 self.info['rate'] = 1 / data_df.time.diff().mean()
163 # Calculate the first and second spectrum
164 meas[field] = self._get_raw_meas(data, **kwargs)
166 if meas[field] == None:
167 meas.pop(field)
170 tmp_df = pd.DataFrame(
171 {'Field': field,
172 'Time': data_df.time.iloc[:max_len],
173 'Vx': data_df.Vx.iloc[:max_len] * kwargs.get('factor', 1e3),
174 'Vy': data_df.Vy.iloc[:max_len] * kwargs.get('factor', 1e3),
175 })
176 all_data = pd.concat([all_data, tmp_df]) #, how='outer')
178 return all_data, meas
180 def _get_raw_meas(self, data, **kwargs):
181 """
182 Returns the RAW Measurement.
183 :param dict data: time signal and info
184 :param bool calculate_second_by_hand: if we calc
185 :return: RAWMeasurement
186 """
187 if len(data['data']) == 0:
188 return None
190 ret = RAW(data,
191 rate=data['info']['rate'],
192 nof_first_spectra=self.info.get('num', 64),
193 calc_first=self.info['first'],
194 calc_second=self.info['second'],
195 downsample=self.info['downsample'],
196 **kwargs.get('raw_kwargs', dict())
197 )
199 if kwargs.get('calculate_second_by_hand', False):
200 ret.calc_second_spectrum(
201 highest_octave=kwargs.get(
202 'highest_octave', ret.avg_spec.freq.max()),
203 minimum_points_in_octave=kwargs.get(
204 'minimum_points_in_octave', 3),
205 peak_filter=kwargs.get(
206 'peak_filter', False),
207 )
209 return ret
211 def cut_fmax(self, fmax):
212 """Cut frequencies higher than maximum allowed frequencies.
214 Args:
215 fmax:
216 """
217 for field, spec in self.measurements.items():
218 # Drop all frequencies larger than fmax
219 s = spec.spectrum.first_spectra_df
220 s.drop(s.query('Frequency > %s' % fmax).index.tolist(),
221 axis=0, inplace=True
222 )
223 freq = spec.spectrum.frequency_span_array
224 freq = freq[freq <= fmax]
225 spec.spectrum.frequency_span_array = freq
227 # Reset First Spectrum Variables
228 spec.spectrum.finalize_first_spectrum(s)
229 spec.avg_spec = pd.DataFrame({
230 'freq': freq,
231 'S': spec.spectrum.first_spectrum_avg.loc[:] # copy
232 })
234 def __equalize_length(self):
235 """
236 Check that data for all fields have same amount of data.
237 :return:
238 """
239 # Check if all fields have same amount of data
240 nofvalues = self.data.groupby('Field').count()
241 if nofvalues.Vx.max() - nofvalues.Vx.min() > 0:
242 fields_with_more_data = nofvalues[
243 nofvalues.Vx > nofvalues.Vx.median()].index.tolist()
244 fields_with_less_data = nofvalues[
245 nofvalues.Vx < nofvalues.Vx.median()].index.tolist()
246 self.logger.warning(
247 'Not all fields have same amount of data (median = %s):\n' % (
248 nofvalues.Vx.quantile(.5))
249 + 'Fields with more data: %s\n' % fields_with_more_data
250 + 'Fields with less data: %s\n' % fields_with_less_data
251 + '%s\n' % nofvalues.Vx.describe()
252 + 'To avoid data loss use the option "equalize_length=False"!')
254 for field in self.data.Field.unique():
255 condition = 'Field == %s' % field
256 if field in fields_with_more_data:
257 # Cut data
258 cut_data = self.data.query(condition).iloc[
259 :int(nofvalues.Vx.median())]
260 self.data.drop(
261 self.data.query(condition).index.tolist(),
262 axis=0, inplace=True)
263 self.data = pd.merge(self.data, cut_data, how='outer')
265 if self.data.query(condition).shape[0] > 0:
266 # recalculate Spectrum
267 rate = 1 / self.data.query(condition).Time.diff().median()
268 self.measurements[field] = self._get_raw_meas(
269 {'data': self.data.query(condition).Vx.to_numpy(),
270 'info': {
271 'Nr': self.info['Nr'],
272 'field': field,
273 'rate': rate,
274 }
275 },
276 )
277 else:
278 # Remove field from measurements
279 self.measurements.pop(field)
280 self.data.drop(
281 self.data.query(condition).index.tolist(),
282 axis=0, inplace=True)
284 self.logger.error(
285 'This field (%s) has more data '
286 '(%s) than median (%s): Cutting data!' % (
287 field,
288 nofvalues.Vx.loc[field],
289 nofvalues.Vx.median()
290 )
291 )
293 if field in fields_with_less_data:
294 self.logger.error(
295 'This field (%s) has less data '
296 '(%s) than median (%s): Deleting field from data!' % (
297 field,
298 nofvalues.Vx.loc[field],
299 nofvalues.Vx.median()
300 )
301 )
302 self.data.drop(
303 self.data.query(condition).index,
304 inplace=True)
305 self.measurements.pop(field)
307 def plot_first_and_second(self, **kwargs):
308 """
309 Args:
310 **kwargs:
311 """
312 sns.set_palette(sns.color_palette("hls", len(self.measurements)))
313 fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1,
314 figsize=(11.69, 8.27), )
316 fields = kwargs.get('fields', self.data.sort_values(
317 'Field')['Field'].unique())
318 for field in fields:
319 # Setting shortcuts
320 try:
321 meas = self.measurements[field]
322 except KeyError:
323 continue
325 s = meas.avg_spec
327 ax1.plot(field, s.S.sum(), 'o')
329 if hasattr(meas.spectrum, 'second_spectra'):
330 marker = "+x*ds"
331 color = "rbgyk"
332 for octave in range(
333 meas.spectrum.number_of_octaves):
334 second = meas.spectrum.second_spectra[octave]
335 ax2.plot(field,
336 second.sum(),
337 "%s%s" % (color[octave], marker[octave])
338 )
340 ax2.set_title('$\\sum$ Second Spectrum')
341 ax2.set_ylabel('$S$ [a.u.]')
343 # Legend for second spectrum
344 field = fields[-1]
345 labels = []
346 if hasattr(self.measurements[field].spectrum, 'second_spectra'):
347 for octave in range(
348 self.measurements[field].spectrum.number_of_octaves):
349 labels.append('%s Octave' % (octave + 1))
350 ax2.legend(labels)
352 ax1.set_title(
353 'm%s: $\\sum$ First Spectrum Noise' % self.info['nr'])
354 ax1.set_ylabel(
355 '$\\langle S_V \\rangle$ [$\\mathrm{V}^2/\\mathrm{Hz}$]')
357 for ax in [ax2, ax1]:
358 ax.set_yscale('log')
359 plt.show()
361 def plot_alpha(self, field, s, ax=None, **kwargs):
362 """
363 Fits a spectrum and plots the slope alpha.
364 :param s: pandas.DataFrame of the spectrum
365 :param ax: axis where to plot the data.
366 :return: linear regression fit
367 """
368 if not ax:
369 fig, ax = plt.subplots()
371 # Fitting alpha
372 s['lnf'] = np.log10(s.freq)
373 s['lnS'] = np.log10(s.S)
375 f = scipy.stats.linregress(
376 s.lnf.iloc[:kwargs.get('fit_range')],
377 s.lnS.iloc[:kwargs.get('fit_range')])
379 plt_obj, = ax.plot(field, -f.slope, kwargs.get('alpha_symbol', 'o'), label='$\\alpha$')
381 return f, plt_obj
383 def plot_info2(self, **kwargs):
384 """
385 Args:
386 **kwargs:
387 """
388 if kwargs.get('fig'):
389 fig = kwargs.get('fig')
390 (ax11, ax12) = kwargs.get('axes')
391 else:
392 fig, (ax11, ax12) = plt.subplots(nrows=2, ncols=1,
393 figsize=(5, 10),)
394 fields = kwargs.get('fields',
395 self.data.sort_values('Field')['Field'].unique())
396 self.plot_noise_chararcteristics(ax11, fields=fields)
397 self.spectra_field_contour(ax12, fields=fields)
399 def plot_noise_chararcteristics(self, ax=None, fields=None, **kwargs):
400 """
401 Args:
402 ax:
403 fields:
404 **kwargs:
405 """
406 if not ax:
407 fig, ax = plt.subplots()
408 fields = kwargs.get('fields', self.data.sort_values('Field')['Field'].unique())
410 par2 = ax.twinx()
412 spectra_df = pd.DataFrame()
413 for i, field in enumerate(fields):
414 # Setting shortcut for Average Spectrum s
415 try:
416 meas = self.measurements[field]
417 if meas == None or len(meas.avg_spec) < 2:
418 continue
419 except KeyError:
420 continue
422 s = meas.avg_spec
424 spectra_df['S%s' % field] = s.S
426 f, alpha = self.plot_alpha(field, s, par2, alpha_symbol="bo",
427 **kwargs)
428 meas.info['alpha'] = -f.slope,
429 meas.info['power'] = f.intercept
431 integral, = ax.plot(field, s.S.sum() * s.freq.diff().iloc[1],
432 'gd', label='Integral')
433 intercept, = ax.plot(field, 10**f.intercept, 'r*',
434 label='$S(f=1~\\mathrm{Hz})$')
435 ax.set_yscale('log')
436 plt.legend(handles=[integral, intercept, alpha], frameon=True,
437 fancybox=True, framealpha=.8)
438 #('Integral', '$S(f=1~\\mathrm{Hz})$', '$\\alpha$'))
440 def plot_info(self, show_field=0, **kwargs):
441 """Plots basic mfn infos
443 Args:
444 show_field:
445 kwargs:
446 """
448 sns.set_palette(sns.color_palette("hls", len(self.measurements)))
449 fig, ((ax11, ax12),
450 (ax21, ax22),
451 (ax31, ax32)) = plt.subplots(nrows=3, ncols=2,
452 figsize=self.style.get('figsize'),
453 gridspec_kw=kwargs.get('gridspec_kw', dict()),
454 # wspace=.3,
455 # hspace=.45)
456 )
457 fields = kwargs.get('fields',
458 self.data.sort_values('Field')['Field'].unique())
459 spectra_df = pd.DataFrame()
460 second = pd.DataFrame()
461 for i, field in enumerate(fields):
462 # Setting shortcut for Average Spectrum s
463 try:
464 meas = self.measurements[field]
465 if meas == None or len(meas.avg_spec) < 2:
466 continue
467 except KeyError:
468 continue
470 s = meas.avg_spec
472 spectra_df['S%s' % field] = s.S
474 ax11.loglog(s.freq, s.S, label='%s T' % field)
476 f, _ = self.plot_alpha(field, s, ax21, **kwargs)
477 meas.info['alpha'] = -f.slope,
478 meas.info['power'] = f.intercept
480 ax22.plot(field, s.S.sum() * s.freq.diff().iloc[1], 'o')
482 if field != float(field):
483 raise ValueError
485 time_df = self.data.query('Field == %f' % field)
486 m = time_df['Vx'].mean()
487 std = time_df['Vx'].std()
488 ax32.errorbar(field, m, yerr=2 * std, elinewidth=2, capsize=4)
490 if hasattr(meas.spectrum, 'second_spectra'):
491 if len(meas.spectrum.second_spectra) > 0:
492 second['S%s' % field] = meas.spectrum.second_spectra[0]
494 # Plot Contour Plot of Averaged Spectra per Field
495 if kwargs.get('plot_field_contour', True):
496 self.spectra_field_contour(ax12, **kwargs)
497 else:
498 field = fields[-1]
499 if kwargs.get('plot_second_sum', False):
500 marker = "+x*ds"
501 color = "rbgyk"
502 for octave in range(
503 self.measurements[
504 field].spectrum.number_of_octaves):
505 ax12.plot(field,
506 self.measurements[
507 field].spectrum.second_spectra[
508 octave].sum(),
509 "%s%s" % (color[octave], marker[octave])
510 )
511 ax12.set_ylabel('$S_2 (f)$')
512 elif kwargs.get('plot_second_freq', False):
513 frequencies = self.measurements[
514 field].spectrum.frequency_span_array_second_spectra
515 smin, smax = second.min().min(), second.max().max()
516 levels = np.logspace(np.log10(smin),
517 np.log10(smax), 20)
518 cs = ax12.contourf(fields,
519 frequencies,
520 second,
521 norm=LogNorm(vmin=smin, vmax=smax),
522 levels=levels,
523 cmap=kwargs.get('contour_cm',
524 self.style.get('default_cm')),
525 )
526 cbar = plt.colorbar(cs, ax=ax12)
527 cbar.set_ticklabels(['%.2e' % _ for _ in levels])
528 ax12.set_ylabel('$f$ [$\\mathrm{Hz}$]')
529 ax12.set_title('Second Spectrum over Field')
530 ax12.set_xlabel('$\\mu_0 H_{ext}$ [$T$]')
531 else:
532 self.data.query("Field == @show_field").plot(
533 x='Time', y='Vx', ax=ax12)
535 # Plot 1/f line in Spectrum
536 xmin, ymin = kwargs.get('xymin', (4e-1, 6e-15))
537 factor, alpha = (1.5, 1)
538 ax11.plot([xmin, xmin * factor],
539 [ymin, ymin / (factor ** alpha)],
540 'k--')
541 ax11.annotate('$1/f$', (xmin * np.sqrt(factor),
542 ymin / (factor ** (alpha / 2))))
544 # Set nice title and axis labels
545 ax11.set_title('Noise PSD')
546 ax11.set_xlabel('$f$ [$\\mathrm{Hz}$]')
547 ax11.set_ylabel('$S_V (f)$ [$\\mathrm{V}^2/\\mathrm{Hz}$]')
549 ax21.set_title('Fit slope')
550 ax21.set_ylabel('$\\alpha$')
552 ax22.set_title('Noise Integral')
554 for ax in [ax22, ax12]:
555 ax.set_yscale('log')
557 self.plot_field_contour(ax31, show_field, **kwargs)
559 ax32.set_title('Hysteresis Loop')
560 ax32.set_ylabel('$V_x$ [$%s\\mathrm{V}$]' % (kwargs.get('unit', 'm')))
561 ax32.set_xlabel('$\\mu_0 H_{ext}$ [$T$]')
563 fig.suptitle('m%d: ' % self.info.get('Nr', 431)
564 + 'Rate=$%d\\,\\mathrm{Hz}$, ' % self.info['rate']
565 + 'Length=$%s\\,\\mathrm{s}$ %s' % (self.info['length'],
566 kwargs.get(
567 'add_info', '')))
569 def plot_field_contour(self, ax31, show_field=0, **kwargs):
570 # Show specific field
571 field_data = self.data.query('Field == %f' % show_field)
572 if field_data.empty:
573 self.logger.error("Can not find timesignal for field (%s)" %
574 show_field)
575 else:
576 try:
577 field_spectra = self.measurements[show_field]
578 if kwargs.get('plot_timesignal', False):
579 field_data.plot(x='Time', y='Vx', legend=False, ax=ax31)
580 ax31.set_title(
581 'Timesignal ($B_{ext} = %.2f\\,\\mathrm{mT}$)' % (
582 show_field * 1e3))
583 ax31.set_ylabel(
584 '$V_x$ [$%s\\mathrm{V}$]' % (kwargs.get('unit', 'm')))
586 # show chunksize of first spectra generated
587 num = field_spectra.spectrum.number_of_first_spectra
588 len_first_spectra = len(field_data) / num
589 for i in range(1, num):
590 x = field_data.iloc[int(len_first_spectra * i)]['Time']
591 ax31.axvline(x, linewidth=.5)
592 elif len(field_spectra.spectrum.first_spectra_df) > 1:
593 s = field_spectra.spectrum.first_spectra_df
594 first_spectra = s.drop('Frequency', axis=1)
595 freq = field_spectra.avg_spec.freq
596 smin = kwargs.get('smin', first_spectra.min().min())
597 smax = kwargs.get('smax', first_spectra.max().max())
598 levels = np.logspace(np.log10(smin),
599 np.log10(smax),
600 kwargs.get('numlevels', 10))
601 cs = ax31.contourf(
602 np.linspace(0, field_spectra.info['Time'],
603 field_spectra.num), # x = Time
604 freq, # y = Frequency
605 first_spectra, # z = all spectra
606 norm=LogNorm(vmin=smin, vmax=smax),
607 levels=levels,
608 cmap=kwargs.get('contour_cm',
609 self.style.get('default_cm')),
610 # levels=kwargs.get('contour_levels'),
611 )
612 cbar = plt.gcf().colorbar(cs, ax=ax31)
613 cbar.set_label('$S_V^t (f)$')
614 cbar.set_ticklabels(['%.2e' % _ for _ in levels])
615 ax31.set_yscale('log')
616 ax31.set_title(
617 'Time resolved PSD ($B_{ext} = '
618 '%.2f\\,\\mathrm{mT}$)' % (show_field * 1e3)
619 )
620 ax31.set_ylabel('$f$ [Hz]')
621 except Exception as e:
622 self.logger.error("Can not plot field (H = %s) data: %s" %
623 (show_field, e))
624 raise ValueError('Error')
625 ax31.set_xlabel('Time [s]')
627 def plot_various_noise_repr(self):
628 sns.set_palette(sns.color_palette("hls", len(self.measurements)))
629 fig, axes = plt.subplots(nrows=4, ncols=4,
630 gridspec_kw=dict(wspace=.3, hspace=.45))
631 for i in range(16):
632 ax = axes[i // 4][i % 4]
633 for field in self.data.sort_values(
634 'Field', ascending=False)['Field'].unique():
635 # Setting shortcut for Average Spectrum s
636 s = self.measurements[field].avg_spec
638 ax.plot(field, s.iloc[i]['S'], 'o')
639 ax.set_title('$S_V (f=%.3f)$' % s.iloc[i].freq)
640 plt.show()
642 def plot_various_noise_fits(self, **kwargs):
643 """
644 Args:
645 **kwargs:
646 """
647 sns.set_palette(sns.color_palette("hls", len(self.measurements)))
648 fig, axes = plt.subplots(nrows=kwargs.get('nrows', 3),
649 ncols=kwargs.get('ncols', 3),
650 gridspec_kw=dict(wspace=.3, hspace=.45),
651 figsize=(11.69, 8.27))
652 for i, fit_range in enumerate(kwargs.get('fit_ranges', range(3, 11))):
653 ax = axes[i // kwargs.get('nrows', 3)][i % kwargs.get('ncols', 3)]
654 for field in self.data.sort_values('Field')['Field'].unique():
655 # Setting shortcut for Average Spectrum s
656 s = self.measurements[field].avg_spec
657 fit = scipy.stats.linregress(s.lnf.iloc[:fit_range],
658 s.lnS.iloc[:fit_range])
660 if kwargs.get('value', 'intercept') == 'intercept':
661 ax.plot(field, 10 ** fit.intercept, 'o')
662 ax.set_yscale('log')
663 else:
664 ax.plot(field, -1 * fit.slope, 'o')
666 ax.set_title('Fit Range = %s' % fit_range)
667 fig.suptitle(kwargs.get('title', 'Noise Fit at 1 Hz'))
669 def plot_compare_timesignals(self, fields, **kwargs):
670 """ Plots a grid of 2x2 subplots to compare timesignals and histogram.
671 Args:
672 fields: list
673 nrows: int Default: len(fields)
674 ncols: int Default: 2
675 sharey: bool Default: True
676 figsize: tuple Default: (11.69, 8.27) in
677 plot_range: int Default: -1
678 """
679 fig, axes = plt.subplots(nrows=kwargs.get('nrows', 1),
680 ncols=kwargs.get('ncols', len(fields)),
681 sharey=kwargs.get('sharey', False),
682 figsize=kwargs.get('figsize', (11.69, 8.27))
683 )
684 twinax = []
685 for i, f in enumerate(fields):
686 ax = axes[i]
687 d = self.data[self.data['Field'] ==
688 f].iloc[:kwargs.get('plot_range', -1)]
689 ax.plot(d['Time'], d['Vx'], **kwargs.get('plot_kw', {}))
690 twinax.append(ax.twiny())
691 sns.distplot(d['Vx'], ax=twinax[i], vertical=True, **kwargs.get('dist_kw', {}))
692 ax.legend(['$V_x (B_{ext}=%.1f\\,\\mathrm{mT})$' % (f * 1e3)])
694 return fig, axes, twinax
696 def plot_various_timesignals(self, fields, **kwargs):
697 """ Plots a grid of 9 subplots with different timesignals
698 Args:
699 fields: list
700 nrows: int Default: len(fields)
701 ncols: int Default: 2
702 sharey: bool Default: True
703 figsize: tuple Default: (11.69, 8.27) in
704 plot_range: int Default: -1
705 """
706 fig, axes = plt.subplots(nrows=kwargs.get('nrows', len(fields)),
707 ncols=kwargs.get('ncols', 2),
708 sharey=kwargs.get('sharey', True),
709 figsize=kwargs.get('figsize', (11.69, 8.27))
710 )
711 for i, f in enumerate(fields):
712 ax = axes[i // kwargs.get('nrows', 3)][i % kwargs.get('ncols', 3)]
713 d = self.data[self.data['Field'] ==
714 f].iloc[:kwargs.get('plot_range', -1)]
715 ax.plot(d['Time'], d['Vx'])
716 ax.legend(['$V_x (B_{ext}=%.1f\\,\\mathrm{mT})$' % (f * 1e3)])
718 def spectra_field_contour(self, ax=None, **kwargs):
719 """
720 Args:
721 ax:
722 **kwargs:
723 """
724 if not ax:
725 fig, ax = plt.subplots()
726 fields = kwargs.get('fields',
727 self.data.sort_values('Field')['Field'].unique())
728 spectra_df = pd.DataFrame()
729 for f in fields:
730 spectra_df[f] = self.measurements[f].avg_spec.S
731 frequencies = self.measurements[f].avg_spec['freq']
732 smin = kwargs.get('smin', spectra_df.min().min())
733 smax = kwargs.get('smax', spectra_df.max().max())
734 levels = np.logspace(np.log10(smin),
735 np.log10(smax),
736 kwargs.get('numlevels', 10))
738 size_fields = len(fields)
739 size_freq = len(frequencies)
740 if spectra_df.shape != (size_freq, size_fields):
741 self.logger.error('Cannot plot contour:\n' + \
742 'spectra_df.shape (%s) != ' + \
743 'size_freq (%s), size_fields (%s)', (
744 spectra_df.shape, size_freq, size_fields))
745 return
747 if size_fields < 2 or size_freq < 2:
748 self.logger.error('Cannot plot contour: ' + \
749 'size_freq (%s), size_fields (%s) incorrect' % (
750 size_freq, size_fields))
751 return
753 cs = ax.contourf(fields * 1e3, frequencies, spectra_df,
754 norm=LogNorm(vmin=smin, vmax=smax),
755 levels=levels,
756 cmap=kwargs.get('contour_cm',
757 self.style.get('default_cm')),
758 )
759 if kwargs.get('cbar', True):
760 cbar = plt.colorbar(cs, ax=ax)
761 cbar.set_ticklabels(['%.2e' % _ for _ in levels])
762 ax.set_yscale('log')
763 ax.set_ylabel('$f$ [$\\mathrm{Hz}$]')
764 ax.set_xlabel('$\\mu_0 H_{ext}$ [$\\mathrm{mT}$]')
765 ax.set_title(kwargs.get('title', 'Noise PSD'))
767 def write(self, file):
768 """Write data to a file.
770 Args:
771 file (str):
772 """
773 if not os.path.exists(file):
774 os.makedirs(file)
776 self.data.to_csv("%s_fields.dat" % file, index=False)
777 for field, m in self.measurements.items():
778 m.write("%s_spectrum_%s" % (file, field))
780 def read(self, file):
781 """Write data to a file.
783 Args:
784 file (str):
785 """
786 self.data = pd.read_csv("%s_fields.dat" % file)
787 for field in self.data['Field'].unique():
788 df = self.data.query('Field == %s' % field)
789 self.measurements[field] = \
790 RAW({'data': df['Vx']},
791 rate=1/df['Time'].diff().median(),
792 )
793 if not self.measurements[field].read("%s_spectrum_%s" % (file,
794 field)):
795 self.measurements.pop(field)
797 def subtract_mean(self, inplace=False, **kwargs):
798 """Subtracts mean value from timesignal.
800 :param inplace: Change variable inside this object?
801 :type inplace: bool
803 :param fields: list of fields to fit
805 :return: normalized timesignal
806 """
808 df_fit = pd.DataFrame()
809 fields = kwargs.get('fields',
810 self.data.sort_values('Field')['Field'].unique())
811 for field in fields:
812 signal = self.data.query('Field == @field').copy()
813 mean = signal.Vx.mean()
814 signal.loc[:, 'normedVx'] = signal.Vx.copy() - mean
816 if inplace:
817 self.data[self.data.Field == field].loc[:, 'normedVx'] = \
818 signal.loc[:, 'normedVx']
820 df_fit = pd.concat([df_fit, signal])
821 return df_fit
823 def plot_multiple_histograms(self, steps=4,
824 fields=pd.Series(dtype=np.float64),
825 factor=1e3,
826 xlim=(-1.3,1.3),
827 **kwargs):
828 """Plots multiple histograms using seaborn ``sns.FacetGrid``
829 with kdeplot.
831 default kwargs for subplots:
832 kde1_kwargs = dict(clip_on=True, shade=True, alpha=1, lw=1.5, bw=.2)
833 kde2_kwargs = kde1_kwargs.update(dict(lw=2, color='w')
834 dist_kwargs = dict(hist=True, norm_hist=True, hist_kws={'alpha': .5})
835 ax_kwargs = dict(y=0, lw=2, clip_on=True).set(xlim=(-1.3,1.3))
836 """
837 if fields.empty:
838 fields = self.data.Field.unique()
840 field_range = abs(fields.max() - fields.min())
841 stepsize = (field_range/steps)
842 for i in range(steps):
843 start = fields.min() + i*stepsize
844 end = start + stepsize
845 df_tmp = self.subtract_mean().query('Field >= @start and Field < @end').copy()
846 df_tmp.loc[:,'normedVx'] *= factor
848 # Initialize the FacetGrid object
849 pal = sns.dark_palette(color="blue", n_colors=10, input='huls')
850 face_kwargs = dict(aspect=5, height=.8, palette=pal)
851 face_kwargs.update(kwargs.get('face_kwargs',{}))
852 kde1_kwargs = dict(clip_on=True, shade=True, alpha=1, lw=1.5, bw=.2)
853 kde1_kwargs.update(kwargs.get('kde1_kwargs',{}))
854 kde2_kwargs = kde1_kwargs.copy()
855 kde2_kwargs.update(dict(lw=2, bw=.2, color='w'))
856 kde2_kwargs.update(kwargs.get('kde2_kwargs',{}))
857 dist_kwargs = dict(hist=True, norm_hist=True, hist_kws={'alpha': .5})
858 dist_kwargs.update(kwargs.get('dist_kwargs',{}))
859 ax_kwargs = dict(y=0, lw=2, clip_on=True)
860 ax_kwargs.update(kwargs.get('ax_kwargs',{}))
862 g = sns.FacetGrid(df_tmp, row="Field", hue="Field", **face_kwargs)
863 g.map(sns.kdeplot, "normedVx", **kde1_kwargs)
864 g.map(sns.kdeplot, "normedVx", **kde2_kwargs)
865 #g.map(sns.distplot, "normedVx", **dist_kwargs)
866 g.map(plt.axhline, **ax_kwargs).set(xlim=xlim)
868 def label(x, color, label):
869 ax = plt.gca()
870 ax.text(0, kwargs.get('adjust_hspace', .3), label, fontweight="bold", color=color,
871 ha="left", va="center", fontsize=16,
872 transform=ax.transAxes)
875 g.map(label, "normedVx")
877 # Set the subplots to overlap
878 g.fig.subplots_adjust(-kwargs.get('adjust_hspace', .3))
880 # Remove axes details that don't play well with overlap
881 g.set_titles("")
882 g.set(yticks=[])#, xticks=[-100 + 50*_ for _ in range(5)])
883 g.despine(bottom=True, left=True)
884 unit = 'mV' if factor == 1e3 else '\\mu{}V' if factor == 1e6 else 'V'
885 plt.gca().set_xlabel('$V_H$ [$\\mathrm{%s}$]' % unit)
887 def load_voltages(self, nr=508, figsize=(16,12)):
888 def load_data(datapath):
889 meas_data = {}
890 meas_info = {}
891 all_data = {}
892 for f in datapath:
893 f_info = self.get_info_from_name(f)
894 sr = f_info['Vin']
895 nr = f_info['nr']
896 meas_info[sr] = f_info
897 meas_data[sr] = pd.read_csv(f, sep='\t')
898 new_df = meas_data[sr]
899 new_df['Vin'] = float(sr[:-2])
900 if nr in all_data.keys():
901 all_data[nr] = pd.concat([all_data[nr], new_df])
902 else:
903 all_data[nr] = new_df
904 return meas_data, meas_info, all_data
906 def calc_PSD(meas_data):
907 meas_obj = {}
908 for sr, data_df in meas_data.items():
909 if len(data_df['Vx']) % 1024:
910 avg = len(data_df['Vx']) // 1024
911 d = data_df['Vx'].iloc[:-(len(data_df['Vx']) % 1024)]
912 else:
913 d = data_df.Vx
915 max_len = len(d)
917 data = {
918 'data': d,
919 'info': {
920 'Nr': meas_info[sr]['nr'],
921 'rate': 1 / data_df.time.diff().mean(),
922 'length': max_len * data_df.time.diff().mean(),
923 }
924 }
926 meas_obj[sr] = ana.RAWMeasurement(data,
927 rate=data['info']['rate'],
928 nof_first_spectra=32,
929 calc_first = True,
930 downsample=False,
931 )
932 return meas_obj
934 def merge_data(meas_obj, cutoff_frequency=.9):
935 diff_voltages = pd.DataFrame()
936 for sr, m in meas_obj.items():
937 s = m.avg_spec
938 s = s[s.freq < cutoff_frequency]
939 if len(s) < 2:
940 continue
941 newdf = pd.DataFrame()
942 newdf['freq'] = s.freq
943 newdf['S'] = s.S
944 newdf['Vin'] = float(sr[:-2])
945 diff_voltages = pd.concat([diff_voltages, newdf])
946 return diff_voltages
948 def plot_PSD_classic(diff_voltages, title, groupby_category='Vin',
949 num=10, style=[['science'], {'context': 'talk', 'style': 'white', 'palette': 'bright',}]):
950 set_style(style)
951 c1 = sns.color_palette("hls", num)
952 sns.set_palette(c1)
953 fig, ax = plt.subplots(figsize=(12,8))
954 #g = sns.relplot(x='freq', y='S', hue='Vin', data=diff_voltages, height=5, kind='line')
955 grouped = diff_voltages.groupby(groupby_category)
956 for group in grouped.groups.keys():
957 grouped.get_group(group).plot(x='freq', y='S', kind='line',
958 loglog=True, ax=ax,
959 label=group,
960 xlabel='Frequency [Hz]',
961 ylabel='$S_{V_H}$ [$\\mathrm{V}^2/\\mathrm{Hz}$]',
962 )
963 ax.set_title(title)
964 #save_plot('m506', 'png')
966 def show_PSD_classic(diff_voltages, title, ax=None, groupby_category='Vin',
967 num=10, style=[['science'], {'context': 'talk', 'style': 'white', 'palette': 'bright',}]):
968 if not ax:
969 fig, ax = plt.subplots(figsize=(12,8))
970 set_style(style)
971 c1 = sns.color_palette("hls", num)
972 sns.set_palette(c1)
973 #g = sns.relplot(x='freq', y='S', hue='Vin', data=diff_voltages, height=5, kind='line')
974 grouped = diff_voltages.groupby(groupby_category)
975 for group in grouped.groups.keys():
976 grouped.get_group(group).plot(x='freq', y='S', kind='line',
977 loglog=True, ax=ax,
978 label=group,
979 xlabel='Frequency [Hz]',
980 ylabel='$S_{V_H}$ [$\\mathrm{V}^2/\\mathrm{Hz}$]',
981 )
982 ax.set_title(title)
983 return ax
986 datapath = glob('./data/MFN/m%s/*' % nr)
987 meas_data, meas_info, all_data = load_data(datapath)
988 meas_obj = calc_PSD(meas_data)
989 diff_voltages = merge_data(meas_obj)
990 fig, ax = plt.subplots(figsize=figsize)
991 show_PSD_classic(diff_voltages, 'm%s: Compare different Voltages' % nr, ax=ax)
993 inset2 = inset_axes(ax, width='100%', height='100%',
994 bbox_to_anchor=(.54, .75, .3, .25),
995 bbox_transform=ax.transAxes)
996 inset3 = inset_axes(ax, width='100%', height='100%',
997 bbox_to_anchor=(.1, .1, .3, .25),
998 bbox_transform=ax.transAxes)
1000 grouped = diff_voltages.groupby('Vin')
1001 for group in grouped.groups.keys():
1002 g = grouped.get_group(group)
1003 fit_area = g.query('freq > %f and freq < %f' % (2e-2, 7e-1))
1004 fit_area['lnf'] = np.log10(fit_area.freq)
1005 fit_area['lnS'] = np.log10(fit_area.S)
1006 fit = scipy.stats.linregress(fit_area.lnf, fit_area.lnS)
1007 intercept, slope = fit.intercept, -fit.slope
1008 voltage = group
1010 inset2.plot(voltage, 10**intercept, 'o')
1011 inset3.plot(voltage, slope, 'o')
1013 inset2.set_xlabel('Voltage [$\mathrm{V}$]')
1014 inset2.set_ylabel('$S_{V_H} (f=1\\;$Hz$)$')
1015 inset2.set_yscale('log')
1017 inset3.set_xlabel('Voltage [$\mathrm{V}$]')
1018 inset3.set_ylabel('$\\alpha$')