diff --git a/config/example-extras.cfg b/config/example-extras.cfg index a99e9dc4..2f2bc4cf 100644 --- a/config/example-extras.cfg +++ b/config/example-extras.cfg @@ -569,6 +569,11 @@ # example, one could set this to "y,x,z" to swap the x and y axes. # It is also possible to negate an axis if the accelerometer # direction is reversed (eg, "x,z,-y"). The default is "x,y,z". +#rate: 3200 +# Output data rate for ADXL345. ADXL345 supports the following data rates: +# 3200, 1600, 800, 400, 200, 100, 50, and 25. Note that it is not recommended +# to change this rate from the default 3200, and rates below 800 will +# considerably affect the quality of resonance measurements. #cs_pin: # The SPI enable pin for the sensor. This parameter must be # provided. @@ -582,6 +587,45 @@ # These optional parameters allow one to customize the SPI settings # used to communicate with the chip. +# Support for resonance testing and automatic input shaper calibration. +# In order to use most of the functionality of this module, additional +# software dependencies must be installed; refer to docs/Measuring_Resonances.md +# for more information. +#[resonance_tester] +#probe_points: +# A list of X,Y,Z coordinates of points (one point per line) to test +# resonances at. At least one point is required. Make sure that all +# points with some safety margin in XY plane (~a few centimeters) are +# reachable by the toolhead. +#accel_chip: +# A name of the accelerometer chip to use for measurements. If adxl345 chip +# was defined without an explicit name, this parameter can simply reference +# it as "accel_chip: adxl345", otherwise an explicit name must be supplied +# as well, e.g. "accel_chip: adxl345 my_chip_name". Either this, or the next +# two parameters must be set. +#accel_chip_x: +#accel_chip_y: +# Names of the accelerometer chips to use for measurements for each of the +# axis. Can be useful, for instance, on bed slinger printer, if two separate +# accelerometers are mounted on the bed (for Y axis) and on the toolhead (for +# X axis). These parameters have the same format as 'accel_chip' parameter. +# Only 'accel_chip' or these two parameters must be provided. +#min_freq: 5 +# Minimum frequency to test for resonances. The default is 5 Hz. +#max_freq: 120 +# Maximum frequency to test for resonances. The default is 120 Hz. +#accel_per_hz: 75 +# This parameter is used to determine which acceleration to use to test a +# specific frequency: accel = accel_per_hz * freq. Higher the value, the +# higher is the energy of the oscillations. Can be set to a lower than the +# default value if the resonances get too strong on the printer. However, +# lower values make measurements of high-frequency resonances less precise. +# The default value is 75 (mm/sec). +#hz_per_sec: 1 +# Determines the speed of the test. When testing all frequencies in range +# [min_freq, max_freq], each second the frequency increases by hz_per_sec. +# Small values make the test slow, and the large values will decrease the +# precision of the test. The default value is 1.0 (Hz/sec == sec^-2). ###################################################################### # Config file helpers diff --git a/docs/G-Codes.md b/docs/G-Codes.md index 8e588a70..b2a35b52 100644 --- a/docs/G-Codes.md +++ b/docs/G-Codes.md @@ -642,14 +642,53 @@ section is enabled: ## Adxl345 Accelerometer Commands -The following command is available when an "adxl345" config section is +The following commands are available when an "adxl345" config section is enabled: - `ACCELEROMETER_MEASURE [CHIP=] [RATE=] [NAME=]`: Starts accelerometer measurements at the requested number of samples per second. If CHIP is not specified it defaults to "default". Valid rates are 25, 50, 100, 200, 400, 800, 1600, - and 3200. If RATE is zero (or not specified) then the current series - of measurements are stopped and the results are written to a file - named `/tmp/adxl345-.csv` where "" is the optional NAME - parameter. If NAME is not specified it defaults to the current time - in "YYYYMMDD_HHMMSS" format. + and 3200. The command works in a start-stop mode: when executed for + the first time, it starts the measurements, next execution stops them. + If RATE is not specified, then the default value is used (either from + `printer.cfg` or `3200` default value). The results of measurements + are written to a file named `/tmp/adxl345-.csv` where "" + is the optional NAME parameter. If NAME is not specified it defaults + to the current time in "YYYYMMDD_HHMMSS" format. +- `ACCELEROMETER_QUERY [CHIP=] [RATE=]`: queries + accelerometer for the current value. If CHIP is not specified it + defaults to "default". If RATE is not specified, the default value is + used. This command is useful to test the connection to the ADXL345 + accelerometer: one of the returned values should be a free-fall + acceleration (+/- some noise of the chip). + +## Resonance Testing Commands + +The following commands are available when a "resonance_tester" config section +is enabled: +- `MEASURE_AXES_NOISE`: Measures and outputs the noise for all axes of all + enabled accelerometer chips. +- `TEST_RESONANCES AXIS= OUTPUT= [NAME=] + [FREQ_START=] [FREQ_END=] [HZ_PER_SEC=]`: + Runs the resonance test in all configured probe points for the requested + axis (X or Y) and measures the acceleration using the accelerometer chips + configured for the respective axis. `OUTPUT` parameter is a comma-separated + list of which outputs will be written. If `raw_data` is requested, then the + raw accelerometer data is written into a file or a series of files + `/tmp/raw_data__[_].csv` with (`_` part of + the name generated only if more than 1 probe point is configured). If + `resonances` is specified, the frequency response is calculated (across + all probe points) and written into `/tmp/resonances__.csv` + file. If unset, OUTPUT defaults to `resonances`, and NAME defaults to + the current time in "YYYYMMDD_HHMMSS" format. +- `SHAPER_CALIBRATE [AXIS=] [NAME=] + [FREQ_START=] [FREQ_END=] [HZ_PER_SEC=]`: + Similarly to `TEST_RESONANCES`, runs the resonance test as configured, and + tries to find the optimal parameters for the input shaper for the requested + axis (or both X and Y axes if `AXIS` parameter is unset). The results of the + tuning are printed to the console, and the frequency responses and the + different input shapers values are written to a CSV file(s) + `/tmp/calibration_data__.csv`. Unless specified, NAME defaults + to the current time in "YYYYMMDD_HHMMSS" format. Note that the suggested + input shaper parameters can be persisted in the config by issuing + `SAVE_CONFIG` command. diff --git a/docs/Measuring_Resonances.md b/docs/Measuring_Resonances.md new file mode 100644 index 00000000..7d0b3263 --- /dev/null +++ b/docs/Measuring_Resonances.md @@ -0,0 +1,348 @@ +Measuring Resonances +==================== + +Klipper has built-in support for ADXL345 accelerometer, which can be used to +measure resonance frequencies of the printer for different axes, and auto-tune +[input shapers](Resonance_Compensation.md) to compensate for resonances. +Note that using ADXL345 requires some soldering and crimping. ADXL345 can be +connected to a Raspberry Pi directly, or to an SPI interface of an MCU +board (it needs to be reasonably fast). + + +Installation instructions +=========================== + +## Wiring + +You need to connect ADXL345 to your Raspberry Pi via SPI. Note that the I2C +connection, which is suggested by ADXL345 documentation, has too low throughput +and **will not work**. The recommended connection scheme: + +| ADXL345 pin | RPi pin | RPi pin name | +|:--:|:--:|:--:| +| 3V3 (or VCC) | 01 | 3.3v DC power | +| GND | 06 | Ground | +| CS | 24 | GPIO08 (SPI0_CE0_N) | +| SDO | 21 | GPIO09 (SPI0_MISO) | +| SDA | 19 | GPIO10 (SPI0_MOSI) | +| SCL | 23 | GPIO11 (SPI0_SCLK) | + +Fritzing wiring diagrams for some of the ADXL345 boards: + +![ADXL345-Rpi](img/adxl345-fritzing.png) + + +Double-check your wiring before powering up the Raspberry Pi to prevent +damaging it or the accelerometer. + +## Mounting the accelerometer + +The accelerometer must be attached to the toolhead. One needs to design a proper +mount that fits their own 3D printer. It is better to align the axes of the +accelerometer with the printer's axes (but if it makes it more convenient, +axes can be swapped - i.e. no need to align X axis with X and so forth - it +should be fine even if Z axis of accelerometer is X axis of the printer, etc.). + +An example of mounting ADXL345 on the SmartEffector: + +![ADXL345 on SmartEffector](img/adxl345-mount.jpg) + +Note that on a bed slinger printer one must design 2 mounts: one for the +toolhead and one for the bed, and run the measurements twice. + +## Software installation + +Note that resonance measurements and shaper auto-calibration require additional +software dependencies not installed by default. You will have to run on your +Raspberry Pi +``` +$ ~/klippy-env/bin/pip install -v numpy +``` +to install `numpy` package. Note that, depending on the performance of the +CPU, it may take *a lot* of time, up to 10-20 minutes. Be patient and wait +for the completion of the installation. On some occasions, if the board has +too little RAM, the installation may fail and you will need to enable swap. + +If installing prerequisites takes too much time or fail for whatever reason, +there is, in principle, another possibility to run a stand-alone script to +automatically tune the input shapers (will be covered later in the guide). + +In order to run stand-alone scripts, one must run the following command to +install the required dependencies (either on Raspberry Pi, or on host, +depending on where the scripts will be executed): +``` +$ sudo apt install python-numpy python-matplotlib +``` + +Afterwards, follow the instructions in the +[RPi Microcontroller document](RPi_microcontroller.md) to setup the +"linux mcu" on the Raspberry Pi. + +Make sure the Linux SPI driver is enabled by running `sudo +raspi-config` and enabling SPI under the "Interfacing options" menu. + +Add the following to the printer.cfg file: +``` +[mcu rpi] +serial: /tmp/klipper_host_mcu + +[adxl345] +cs_pin: rpi:None + +[resonance_tester] +accel_chip: adxl345 +probe_points: + 100,100,20 # an example +``` +It is advised to start with 1 probe point, in the middle of the print bed, +slightly above it. + +Restart Klipper via the `RESTART` command. + +Measuring the resonances +=========================== + +## Checking the setup + +Now you can test a connection. In Octoprint, run `ACCELEROMETER_QUERY`. You +should see the current measurements from the accelerometer, including the +free-fall acceleration, e.g. +``` +Recv: // adxl345 values (x, y, z): 470.719200, 941.438400, 9728.196800 +``` + +Try running `MEASURE_AXES_NOISE` in Octoprint, you should get some baseline +numbers for the noise of accelerometer on the axes (should be somewhere +in the range of ~1-100). Note that this feature will not be available if +`numpy` package was not installed (see +[Software installation](#software-installation) for more details). + +## Measuring the resonances + +Now you can run some real-life tests. In `printer.cfg` add or replace the +following values: +``` +[printer] +max_accel: 7000 +max_accel_to_decel: 7000 +``` +(after you are done with the measurements, revert these values to their old, +or the newly suggested values). Also, if you have enabled input shaper already, +you will need to disable it prior to this test as follows: +``` +SET_INPUT_SHAPER SHAPER_FREQ_X=0 SHAPER_FREQ_Y=0 +``` +as it is not valid to run the resonance testing with the input shaper enabled. + +Run the following command: +``` +TEST_RESONANCES AXIS=X +``` +Note that it will create vibrations on X axis. If that works, run for Y axis +as well: +``` +TEST_RESONANCES AXIS=Y +``` +This will generate 2 CSV files (`/tmp/resonances_x_*.csv` and +`/tmp/resonances_y_*.csv`). + +Note that the commands above require `numpy` to be installed installed. If you +haven't installed it, you can instead pass `OUTPUT=raw_data` argument to the +above commands (2 files `/tmp/raw_data_x_*.csv` and `/tmp/raw_data_y_*.csv` +will be written). One can then run stand-alone scripts on Raspberry Pi +(specify the correct file name on the command line): +``` +$ ~/klipper/scripts/graph_accelerometer.py /tmp/raw_data_x_*.csv -o /tmp/resonances_x.png +$ ~/klipper/scripts/calibrate_shaper.py /tmp/raw_data_x_*.csv -o /tmp/shaper_calibrate_x.png +``` +or copy the data to the host and run the scripts there. See +[Offline processing of the accelerometer data](#offline-processing-of-the-accelerometer-data) +section for more details. + +**Attention!** Be sure to observe the printer for the first time, to make sure +the vibrations do not become too violent (`M112` command can be used to abort +the test in case of emergency; hopefully it will not come to this though). +If the vibrations do get too strong, you can attempt to specify a lower than the +default value for `accel_per_hz` parameter in `[resonance_tester]` section, e.g. +``` +[resonance_tester] +accel_chip: adxl345 +accel_per_hz: 50 # default is 75 +probe_points: ... +``` + +Generated CSV files show power spectral density of the vibrations depending on the +frequency. Usually, the charts generated from these CSV files are relatively easy +to read, with the peaks corresponding to the resonance frequencies: + +![Resonances](img/test-resonances-x.png) + +The chart above shows the resonances for X axis at approx. 50 Hz, 56 Hz, 63 Hz, +80 Hz and 104 Hz and one cross-resonance for Y axis at ~ 56 Hz. From this, one +can derive that a good input shaper config in this case could be `2hump_ei` at +around `shaper_freq_y = 45` (Hz): + +|![2-hump EI shaper](img/2hump_ei_65hz.png)| +|:--:| +|Input Shaper response to vibrations, lower is better.| + +Note that the smaller resonance at 104 Hz requires less of vibration suppression +(if at all). + +## Input Shaper auto-calibration + +Besides manually choosing the appropriate parameters for the input shaper +feature, it is also possible to run an experimental auto-tuning for the +input shaper. + +In order to attempt to measure the resonance frequencies and automatically +determine the best parameters for `[input_shaper]`, run the following command +via Octoprint terminal: +``` +SHAPER_CALIBRATE +``` + +This will test all frequencies in range 5 Hz - 120 Hz and generate +the csv output (`/tmp/calibration_data_*.csv` by default) for the frequency +response and the suggested input shapers. You will also get the suggested +frequencies for each input shaper, as well as which input shaper is recommended +for your setup, on Octoprint console. For example: + +![Resonances](img/calibrate-y.png) +``` +Fitted shaper 'zv' frequency = 56.7 Hz (vibrations = 23.2%) +Fitted shaper 'mzv' frequency = 52.9 Hz (vibrations = 10.9%) +Fitted shaper 'ei' frequency = 62.0 Hz (vibrations = 8.9%) +Fitted shaper '2hump_ei' frequency = 59.0 Hz (vibrations = 4.9%) +Fitted shaper '3hump_ei' frequency = 65.0 Hz (vibrations = 3.3%) +Recommended shaper_type_y = 2hump_ei, shaper_freq_y = 59.0 Hz +``` +If you agree with the suggested parameters, you can execute `SAVE_CONFIG` +now to save them and restart the Klipper. + + +If your printer is a bed slinger printer, you will need to repeat the +measurements twice: measure the resonances of X axis with the accelerometer +attached to the toolhead and the resonances of Y axis - to the bed (the usual +bed slinger setup). In this case, you can specify the axis you want to run the +test for (by default the test is performed for both axes): +``` +SHAPER_CALIBRATE AXIS=Y +``` + +You can execute `SAVE_CONFIG` twice - after calibrating each axis. + +However, you can connect two accelerometers simultaneously, though they must be +connected to different boards (say, to an RPi and printer MCU board), or to two +different physical SPI interfaces on the same board (rarely available). +Then they can be configured in the following manner: +``` +[adxl345 adxl345_x] +# Assuming adxl345_x is connected to an RPi +cs_pin: rpi:None + +[adxl345 adxl345_y] +# Assuming adxl345_y is connected to a printer MCU board +cs_pin: ... # Printer board SPI chip select (CS) pin + +[resonance_tester] +accel_chip_x: adxl345_x +accel_chip_y: adxl345_y +probe_points: ... +``` +then one can simply run `SHAPER_CALIBRATE` without specifying an axis to +calibrate the input shaper for both axes in one go. + +After the autocalibration is finished, you will still need to choose the +`max_accel` value that does not create too much smoothing in the printed +parts. Follow [this](Resonance_Compensation.md#selecting-max_accel) part of +the input shaper tuning guide and print the test model. + +## Input Shaper re-calibration + +`SHAPER_CALIBRATE` command can be also used to re-calibrate the input shaper in +the future, especially if some changes to the printer that can affect its +kinematics are made. One can either re-run the full calibration using +`SHAPER_CALIBRATE` command, or restrict the auto-calibration to a single axis by +supplying `AXIS=` parameter, like +``` +SHAPER_CALIBRATE AXIS=X +``` + +**Warning!** It is not advisable to run the shaper autocalibration very +frequently (e.g. before every print, or every day). In order to determine +resonance frequencies, autocalibration creates intensive vibrations on each of +the axes. Generally, 3D printers are not designed to withstand a prolonged +exposure to vibrations near the resonance frequencies. Doing so may increase +wear of the printer components and reduce their lifespan. There is also an +increased risk of some parts unscrewing or becoming loose. Always check that +all parts of the printer (including the ones that may normally not move) are +securely fixed in place after each auto-tuning. + +Also, due to some noise in measurements, it is possible the the tuning results +will be slightly different from one calibration run to another one. Still, it +is not expected that the resulting print quality will be affected too much. +However, it is still advised to double-check the suggested parameters, and +print some test prints before using them to confirm they are good. + +## Offline processing of the accelerometer data + +It is possible to generate the raw accelerometer data and process it offline +(e.g. on a host machine), for example to find resonances. In order to do so, +run the following command via Octoprint terminal: +``` +TEST_RESONANCES AXIS=X OUTPUT=raw_data +``` +(specify the desired test axis and the desired template for the raw +accelerometer output, the data will be written into `/tmp` directory). + +The raw data can also be obtained by running the command `ACCELEROMETER_MEASURE` +command twice during some normal printer activity - first to start the +measurements, and then to stop them and write the output file. Refer to +[G-Codes](G-Codes.md#adxl345-accelerometer-commands) for more details. + +The data can be processed later by the following scripts: +`scripts/graph_accelerometer.py` and `scripts/calibrate_shaper.py`. Both +of them accept one or several raw csv files as the input depending on the +mode. The graph_accelerometer.py script supports several modes of operation: + * plotting raw accelerometer data (use `-r` parameter), only 1 input is + supported; + * plotting a frequency response (no extra parameters required), if multiple + inputs are specified, the average frequency response is computed; + * comparison of the frequency response between several inputs (use `-c` + parameter); you can additionally specify which accelerometer axis to + consider via `-a x`, `-a y` or `-a z` parameter (if none specified, + the sum of vibrations for all axes is used); + * plotting the spectrogram (use `-s` parameter), only 1 input is supported; + you can additionally specify which accelerometer axis to consider via + `-a x`, `-a y` or `-a z` parameter (if none specified, the sum of vibrations + for all axes is used). + +For example, +``` +$ ~/klipper/scripts/graph_accelerometer.py /tmp/raw_data_x_*.csv -o /tmp/resonances_x.png -c -a z +``` +will plot the comparison of several `/tmp/raw_data_x_*.csv` files for Z axis to +`/tmp/resonances_x.png` file. + +The shaper_calibrate.py script accepts 1 or several inputs and can run automatic +tuning of the input shaper and suggest the best parameters that work well for +all provided inputs. It prints the suggested parameters to the console, and can +additionally generate the chart if `-o output.png` parameter is provided, or +the CSV file if `-c output.csv` parameter is specified. + +Providing several inputs to shaper_calibrate.py script can be useful if running +some advanced tuning of the input shapers, for example: + * Running `TEST_RESONANCES AXIS=X OUTPUT=raw_data` (and `Y` axis) for a single + axis twice on a bed slinger printer with the accelerometer attached to the + toolhead the first time, and the accelerometer attached to the bed the + second time in order to detect axes cross-resonances and attempt to cancel + them with input shapers. + * Running `TEST_RESONANCES AXIS=Y OUTPUT=raw_data` twice on a bed slinger with + a glass bed and a magnetic surfaces (which is lighter) to find the input + shaper parameters that work well for any print surface configuration. + * Combining the resonance data from multiple test points. + * Combining the resonance data from 2 axis (e.g. on a bed slinger printer + to configure X-axis input_shaper from both X and Y axes resonances to + cancel vibrations of the *bed* in case the nozzle 'catches' a print when + moving in X axis direction). diff --git a/docs/Resonance_Compensation.md b/docs/Resonance_Compensation.md index 4d5fd3aa..c2706549 100644 --- a/docs/Resonance_Compensation.md +++ b/docs/Resonance_Compensation.md @@ -195,7 +195,10 @@ A few notes on shaper selection: ## Selecting max_accel -You should have a printed test for the shaper you chose from the previous step. +You should have a printed test for the shaper you chose from the previous step +(if you don't, print the test model with the pressure advance disabled +`SET_PRESSURE_ADVANCE ADVANCE=0` and with the tuning tower enabled as +`TUNING_TOWER COMMAND=SET_VELOCITY_LIMIT PARAMETER=ACCEL START=1250 FACTOR=100 BAND=5`). Note that at very high accelerations, depending on the resonance frequency and the input shaper you chose (e.g. EI shaper creates more smoothing than MZV), input shaping may cause too much smoothing and rounding of the parts. So, diff --git a/docs/img/2hump_ei_65hz.png b/docs/img/2hump_ei_65hz.png new file mode 100644 index 00000000..56ad2f79 Binary files /dev/null and b/docs/img/2hump_ei_65hz.png differ diff --git a/docs/img/adxl345-fritzing.png b/docs/img/adxl345-fritzing.png new file mode 100644 index 00000000..018ba5d6 Binary files /dev/null and b/docs/img/adxl345-fritzing.png differ diff --git a/docs/img/adxl345-mount.jpg b/docs/img/adxl345-mount.jpg new file mode 100644 index 00000000..cc05d885 Binary files /dev/null and b/docs/img/adxl345-mount.jpg differ diff --git a/docs/img/calibrate-y.png b/docs/img/calibrate-y.png new file mode 100644 index 00000000..520eecc1 Binary files /dev/null and b/docs/img/calibrate-y.png differ diff --git a/docs/img/test-resonances-x.png b/docs/img/test-resonances-x.png new file mode 100644 index 00000000..e69336c6 Binary files /dev/null and b/docs/img/test-resonances-x.png differ diff --git a/klippy/extras/adxl345.py b/klippy/extras/adxl345.py index fb2928c7..c29d5451 100644 --- a/klippy/extras/adxl345.py +++ b/klippy/extras/adxl345.py @@ -3,7 +3,7 @@ # Copyright (C) 2020 Kevin O'Connor # # This file may be distributed under the terms of the GNU GPLv3 license. -import logging, time, collections +import logging, time, collections, multiprocessing, os from . import bus # ADXL345 registers @@ -28,11 +28,10 @@ Accel_Measurement = collections.namedtuple( # Sample results class ADXL345Results: def __init__(self): + self.raw_samples = None self.samples = [] self.drops = self.overflows = 0 self.time_per_sample = self.start_range = self.end_range = 0. - def get_samples(self): - return self.samples def get_stats(self): return ("drops=%d,overflows=%d" ",time_per_sample=%.9f,start_range=%.6f,end_range=%.6f" @@ -42,31 +41,57 @@ class ADXL345Results: start1_time, start2_time, end1_time, end2_time): if not raw_samples or not end_sequence: return - (x_pos, x_scale), (y_pos, y_scale), (z_pos, z_scale) = axes_map + self.axes_map = axes_map + self.raw_samples = raw_samples self.overflows = overflows + self.start2_time = start2_time self.start_range = start2_time - start1_time self.end_range = end2_time - end1_time - total_count = (end_sequence - 1) * 8 + len(raw_samples[-1][1]) // 6 + self.total_count = (end_sequence - 1) * 8 + len(raw_samples[-1][1]) // 6 total_time = end2_time - start2_time - self.time_per_sample = time_per_sample = total_time / total_count - seq_to_time = time_per_sample * 8. - self.samples = samples = [None] * total_count + self.time_per_sample = time_per_sample = total_time / self.total_count + self.seq_to_time = time_per_sample * 8. + actual_count = sum([len(data)//6 for _, data in raw_samples]) + self.drops = self.total_count - actual_count + def decode_samples(self): + if not self.raw_samples: + return self.samples + (x_pos, x_scale), (y_pos, y_scale), (z_pos, z_scale) = self.axes_map actual_count = 0 - for seq, data in raw_samples: + self.samples = samples = [None] * self.total_count + for seq, data in self.raw_samples: d = bytearray(data) count = len(data) sdata = [(d[i] | (d[i+1] << 8)) - ((d[i+1] & 0x80) << 9) for i in range(0, count-1, 2)] - seq_time = start2_time + seq * seq_to_time + seq_time = self.start2_time + seq * self.seq_to_time for i in range(count//6): - samp_time = seq_time + i * time_per_sample + samp_time = seq_time + i * self.time_per_sample x = sdata[i*3 + x_pos] * x_scale y = sdata[i*3 + y_pos] * y_scale z = sdata[i*3 + z_pos] * z_scale samples[actual_count] = Accel_Measurement(samp_time, x, y, z) actual_count += 1 del samples[actual_count:] - self.drops = total_count - actual_count + return self.samples + def write_to_file(self, filename): + def write_impl(): + try: + # Try to re-nice writing process + os.nice(20) + except: + pass + f = open(filename, "w") + f.write("##%s\n#time,accel_x,accel_y,accel_z\n" % ( + self.get_stats(),)) + samples = self.samples or self.decode_samples() + for t, accel_x, accel_y, accel_z in samples: + f.write("%.6f,%.6f,%.6f,%.6f\n" % ( + t, accel_x, accel_y, accel_z)) + f.close() + write_proc = multiprocessing.Process(target=write_impl) + write_proc.daemon = True + write_proc.start() # Printer class that controls measurments class ADXL345: @@ -80,6 +105,9 @@ class ADXL345: if len(axes_map) != 3 or any([a.strip() not in am for a in axes_map]): raise config.error("Invalid adxl345 axes_map parameter") self.axes_map = [am[a.strip()] for a in axes_map] + self.data_rate = config.getint('rate', 3200) + if self.data_rate not in QUERY_RATES: + raise config.error("Invalid rate parameter: %d" % (self.data_rate,)) # Measurement storage (accessed from background thread) self.raw_samples = [] self.last_sequence = 0 @@ -104,9 +132,14 @@ class ADXL345: gcode.register_mux_command("ACCELEROMETER_MEASURE", "CHIP", name, self.cmd_ACCELEROMETER_MEASURE, desc=self.cmd_ACCELEROMETER_MEASURE_help) + gcode.register_mux_command("ACCELEROMETER_QUERY", "CHIP", name, + self.cmd_ACCELEROMETER_QUERY, + desc=self.cmd_ACCELEROMETER_QUERY_help) if name == "default": gcode.register_mux_command("ACCELEROMETER_MEASURE", "CHIP", None, self.cmd_ACCELEROMETER_MEASURE) + gcode.register_mux_command("ACCELEROMETER_QUERY", "CHIP", None, + self.cmd_ACCELEROMETER_QUERY) def _build_config(self): self.query_adxl345_cmd = self.mcu.lookup_command( "query_adxl345 oid=%c clock=%u rest_ticks=%u", @@ -128,7 +161,7 @@ class ADXL345: sequence += 0x10000 self.last_sequence = sequence raw_samples = self.raw_samples - if len(raw_samples) >= 200000: + if len(raw_samples) >= 300000: # Avoid filling up memory with too many samples return raw_samples.append((sequence, params['data'])) @@ -138,8 +171,7 @@ class ADXL345: sequence += 0x10000 return sequence def start_measurements(self, rate=None): - if rate is None: - rate = 3200 + rate = rate or self.data_rate # Verify chip connectivity params = self.spi.spi_transfer([REG_DEVID | REG_MOD_READ, 0x00]) response = bytearray(params['response']) @@ -190,33 +222,47 @@ class ADXL345: self.samples_start1, self.samples_start2, end1_time, end2_time) logging.info("ADXL345 finished %d measurements: %s", - len(res.get_samples()), res.get_stats()) + res.total_count, res.get_stats()) return res def end_query(self, name): if not self.query_rate: return res = self.finish_measurements() # Write data to file - f = open("/tmp/adxl345-%s.csv" % (name,), "w") - f.write("##%s\n#time,accel_x,accel_y,accel_z\n" % (res.get_stats(),)) - for t, accel_x, accel_y, accel_z in res.get_samples(): - f.write("%.6f,%.6f,%.6f,%.6f\n" % (t, accel_x, accel_y, accel_z)) - f.close() + filename = "/tmp/adxl345-%s.csv" % (name,) + res.write_to_file(filename) cmd_ACCELEROMETER_MEASURE_help = "Start/stop accelerometer" def cmd_ACCELEROMETER_MEASURE(self, gcmd): - rate = gcmd.get_int("RATE", 0) - if not rate: + if self.query_rate: name = gcmd.get("NAME", time.strftime("%Y%m%d_%H%M%S")) if not name.replace('-', '').replace('_', '').isalnum(): raise gcmd.error("Invalid adxl345 NAME parameter") self.end_query(name) gcmd.respond_info("adxl345 measurements stopped") - elif self.query_rate: - raise gcmd.error("adxl345 already running") - elif rate not in QUERY_RATES: - raise gcmd.error("Not a valid adxl345 query rate") else: + rate = gcmd.get_int("RATE", self.data_rate) + if rate not in QUERY_RATES: + raise gcmd.error("Not a valid adxl345 query rate: %d" % (rate,)) self.start_measurements(rate) + gcmd.respond_info("adxl345 measurements started") + cmd_ACCELEROMETER_QUERY_help = "Query accelerometer for the current values" + def cmd_ACCELEROMETER_QUERY(self, gcmd): + if self.query_rate: + raise gcmd.error("adxl345 measurements in progress") + self.start_measurements() + reactor = self.printer.get_reactor() + eventtime = starttime = reactor.monotonic() + while not self.raw_samples: + eventtime = reactor.pause(eventtime + .1) + if eventtime > starttime + 3.: + # Try to shutdown the measurements + self.finish_measurements() + raise gcmd.error("Timeout reading adxl345 data") + result = self.finish_measurements() + values = result.decode_samples() + _, accel_x, accel_y, accel_z = values[-1] + gcmd.respond_info("adxl345 values (x, y, z): %.6f, %.6f, %.6f" % ( + accel_x, accel_y, accel_z)) def load_config(config): return ADXL345(config) diff --git a/klippy/extras/input_shaper.py b/klippy/extras/input_shaper.py index 35e7ff63..ceca89f8 100644 --- a/klippy/extras/input_shaper.py +++ b/klippy/extras/input_shaper.py @@ -25,14 +25,12 @@ class InputShaper: , 'ei': ffi_lib.INPUT_SHAPER_EI , '2hump_ei': ffi_lib.INPUT_SHAPER_2HUMP_EI , '3hump_ei': ffi_lib.INPUT_SHAPER_3HUMP_EI} - shaper_type = config.getchoice('shaper_type', self.shapers, None) - if shaper_type is None: - self.shaper_type_x = config.getchoice( - 'shaper_type_x', self.shapers, 'mzv') - self.shaper_type_y = config.getchoice( - 'shaper_type_y', self.shapers, 'mzv') - else: - self.shaper_type_x = self.shaper_type_y = shaper_type + shaper_type = config.get('shaper_type', 'mzv') + self.shaper_type_x = config.getchoice( + 'shaper_type_x', self.shapers, shaper_type) + self.shaper_type_y = config.getchoice( + 'shaper_type_y', self.shapers, shaper_type) + self.saved_shaper_freq_x = self.saved_shaper_freq_y = 0. self.stepper_kinematics = [] self.orig_stepper_kinematics = [] # Register gcode commands @@ -86,6 +84,24 @@ class InputShaper: , shaper_type_x, shaper_type_y , shaper_freq_x, shaper_freq_y , damping_ratio_x, damping_ratio_y) + def disable_shaping(self): + if (self.saved_shaper_freq_x or self.saved_shaper_freq_y) and not ( + self.shaper_freq_x or self.shaper_freq_y): + # Input shaper is already disabled + return + self.saved_shaper_freq_x = self.shaper_freq_x + self.saved_shaper_freq_y = self.shaper_freq_y + self._set_input_shaper(self.shaper_type_x, self.shaper_type_y, 0., 0., + self.damping_ratio_x, self.damping_ratio_y) + def enable_shaping(self): + saved = self.saved_shaper_freq_x or self.saved_shaper_freq_y + if saved: + self._set_input_shaper(self.shaper_type_x, self.shaper_type_y, + self.saved_shaper_freq_x, + self.saved_shaper_freq_y, + self.damping_ratio_x, self.damping_ratio_y) + self.saved_shaper_freq_x = self.saved_shaper_freq_y = 0. + return saved cmd_SET_INPUT_SHAPER_help = "Set cartesian parameters for input shaper" def cmd_SET_INPUT_SHAPER(self, gcmd): damping_ratio_x = gcmd.get_float( diff --git a/klippy/extras/resonance_tester.py b/klippy/extras/resonance_tester.py new file mode 100644 index 00000000..70bea4a1 --- /dev/null +++ b/klippy/extras/resonance_tester.py @@ -0,0 +1,303 @@ +# A utility class to test resonances of the printer +# +# Copyright (C) 2020 Dmitry Butyugin +# +# This file may be distributed under the terms of the GNU GPLv3 license. +import logging, math, os, time +from . import shaper_calibrate + +def _parse_probe_points(config): + points = config.get('probe_points').split('\n') + try: + points = [line.split(',', 2) for line in points if line.strip()] + return [[float(coord.strip()) for coord in p] for p in points] + except: + raise config.error("Unable to parse probe_points in %s" % ( + config.get_name())) + +class VibrationPulseTest: + def __init__(self, config): + printer = config.get_printer() + self.gcode = printer.lookup_object('gcode') + self.min_freq = config.getfloat('min_freq', 5., minval=1.) + self.max_freq = config.getfloat('max_freq', 120., + minval=self.min_freq, maxval=200.) + self.accel_per_hz = config.getfloat('accel_per_hz', 75.0, above=0.) + self.hz_per_sec = config.getfloat('hz_per_sec', 1., + minval=0.1, maxval=2.) + + self.probe_points = _parse_probe_points(config) + def get_supported_axes(self): + return ['x', 'y'] + def get_start_test_points(self): + return self.probe_points + def prepare_test(self, toolhead, gcmd): + self.freq_start = gcmd.get_float("FREQ_START", self.min_freq, minval=1.) + self.freq_end = gcmd.get_float("FREQ_END", self.max_freq, + minval=self.freq_start, maxval=200.) + self.hz_per_sec = gcmd.get_float("HZ_PER_SEC", self.hz_per_sec, + above=0., maxval=2.) + # Attempt to adjust maximum acceleration and acceleration to + # deceleration based on the maximum test frequency. + max_accel = self.freq_end * self.accel_per_hz + toolhead.cmd_SET_VELOCITY_LIMIT(self.gcode.create_gcode_command( + "SET_VELOCITY_LIMIT", "SET_VELOCITY_LIMIT", + {"ACCEL": max_accel, "ACCEL_TO_DECEL": max_accel})) + def run_test(self, toolhead, axis, gcmd): + X, Y, Z, E = toolhead.get_position() + if axis not in self.get_supported_axes(): + raise gcmd.error("Test axis '%s' is not supported", axis) + vib_dir = (1, 0) if axis == 'x' else (0., 1.) + sign = 1. + freq = self.freq_start + gcmd.respond_info("Testing frequency %.0f Hz" % (freq,)) + _, max_accel = toolhead.get_max_velocity() + while freq <= self.freq_end + 0.000001: + t_seg = .25 / freq + accel = min(self.accel_per_hz * freq, max_accel) + V = accel * t_seg + toolhead.cmd_M204(self.gcode.create_gcode_command( + "M204", "M204", {"S": accel})) + L = .5 * accel * t_seg**2 + nX = X + sign * vib_dir[0] * L + nY = Y + sign * vib_dir[1] * L + toolhead.move([nX, nY, Z, E], V) + toolhead.move([X, Y, Z, E], V) + sign = -sign + old_freq = freq + freq += 2. * t_seg * self.hz_per_sec + if math.floor(freq) > math.floor(old_freq): + gcmd.respond_info("Testing frequency %.0f Hz" % (freq,)) + +class ResonanceTester: + def __init__(self, config): + self.printer = config.get_printer() + self.move_speed = config.getfloat('move_speed', 50., above=0.) + self.test = VibrationPulseTest(config) + if not config.get('accel_chip_x', None): + self.accel_chip_names = [('xy', config.get('accel_chip').strip())] + else: + self.accel_chip_names = [ + ('x', config.get('accel_chip_x').strip()), + ('y', config.get('accel_chip_y').strip())] + if self.accel_chip_names[0][1] == self.accel_chip_names[1][1]: + self.accel_chip_names = [('xy', self.accel_chip_names[0][1])] + + self.gcode = self.printer.lookup_object('gcode') + self.gcode.register_command("MEASURE_AXES_NOISE", + self.cmd_MEASURE_AXES_NOISE) + self.gcode.register_command("TEST_RESONANCES", + self.cmd_TEST_RESONANCES) + self.gcode.register_command("SHAPER_CALIBRATE", + self.cmd_SHAPER_CALIBRATE) + self.printer.register_event_handler("klippy:connect", self.connect) + + def connect(self): + self.accel_chips = [ + (axis, self.printer.lookup_object(chip_name)) + for axis, chip_name in self.accel_chip_names] + + def cmd_TEST_RESONANCES(self, gcmd): + toolhead = self.printer.lookup_object('toolhead') + # Parse parameters + self.test.prepare_test(toolhead, gcmd) + if len(self.test.get_supported_axes()) > 1: + axis = gcmd.get("AXIS").lower() + else: + axis = gcmd.get("AXIS", self.test.get_supported_axes()[0]).lower() + if axis not in self.test.get_supported_axes(): + raise gcmd.error("Unsupported axis '%s'" % (axis,)) + + outputs = gcmd.get("OUTPUT", "resonances").lower().split(',') + for output in outputs: + if output not in ['resonances', 'raw_data']: + raise gcmd.error("Unsupported output '%s', only 'resonances'" + " and 'raw_data' are supported" % (output,)) + if not outputs: + raise gcmd.error("No output specified, at least one of 'resonances'" + " or 'raw_data' must be set in OUTPUT parameter") + name_suffix = gcmd.get("NAME", time.strftime("%Y%m%d_%H%M%S")) + if not self.is_valid_name_suffix(name_suffix): + raise gcmd.error("Invalid NAME parameter") + csv_output = 'resonances' in outputs + raw_output = 'raw_data' in outputs + + # Setup calculation of resonances + if csv_output: + helper = shaper_calibrate.ShaperCalibrate(self.printer) + + currentPos = toolhead.get_position() + Z = currentPos[2] + E = currentPos[3] + + calibration_points = self.test.get_start_test_points() + data = None + for point in calibration_points: + toolhead.manual_move(point, self.move_speed) + if len(calibration_points) > 1: + gcmd.respond_info( + "Probing point (%.3f, %.3f, %.3f)" % tuple(point)) + toolhead.wait_moves() + toolhead.dwell(0.500) + gcmd.respond_info("Testing axis %s" % axis.upper()) + + for chip_axis, chip in self.accel_chips: + if axis in chip_axis or chip_axis in axis: + chip.start_measurements() + # Generate moves + self.test.run_test(toolhead, axis, gcmd) + raw_values = [] + for chip_axis, chip in self.accel_chips: + if axis in chip_axis or chip_axis in axis: + results = chip.finish_measurements() + if raw_output: + raw_name = self.get_filename( + 'raw_data', name_suffix, axis, + point if len(calibration_points) > 1 else None) + results.write_to_file(raw_name) + gcmd.respond_info( + "Writing raw accelerometer data to %s file" % ( + raw_name,)) + raw_values.append((chip_axis, results)) + if not csv_output: + continue + for chip_axis, chip_values in raw_values: + gcmd.respond_info("%s-axis accelerometer stats: %s" % ( + chip_axis, chip_values.get_stats(),)) + if not chip_values: + raise gcmd.error( + "%s-axis accelerometer measured no data" % ( + chip_axis,)) + new_data = helper.process_accelerometer_data(chip_values) + data = data.join(new_data) if data else new_data + if csv_output: + csv_name = self.save_calibration_data('resonances', name_suffix, + helper, axis, data) + gcmd.respond_info( + "Resonances data written to %s file" % (csv_name,)) + + def cmd_SHAPER_CALIBRATE(self, gcmd): + toolhead = self.printer.lookup_object('toolhead') + # Parse parameters + self.test.prepare_test(toolhead, gcmd) + axis = gcmd.get("AXIS", None) + if not axis: + calibrate_axes = self.test.get_supported_axes() + elif axis.lower() not in self.test.get_supported_axes(): + raise gcmd.error("Unsupported axis '%s'" % (axis,)) + else: + calibrate_axes = [axis.lower()] + + name_suffix = gcmd.get("NAME", time.strftime("%Y%m%d_%H%M%S")) + if not self.is_valid_name_suffix(name_suffix): + raise gcmd.error("Invalid NAME parameter") + + # Setup shaper calibration + helper = shaper_calibrate.ShaperCalibrate(self.printer) + + input_shaper = self.printer.lookup_object('input_shaper', None) + if input_shaper is not None: + input_shaper.disable_shaping() + gcmd.respond_info("Disabled [input_shaper] for calibration") + + currentPos = toolhead.get_position() + Z = currentPos[2] + E = currentPos[3] + calibration_data = {axis: None for axis in calibrate_axes} + + calibration_points = self.test.get_start_test_points() + for point in calibration_points: + toolhead.manual_move(point, self.move_speed) + if len(calibration_points) > 1: + gcmd.respond_info( + "Probing point (%.3f, %.3f, %.3f)" % tuple(point)) + for axis in calibrate_axes: + toolhead.wait_moves() + toolhead.dwell(0.500) + gcmd.respond_info("Testing axis %s" % axis.upper()) + + for chip_axis, chip in self.accel_chips: + if axis in chip_axis or chip_axis in axis: + chip.start_measurements() + # Generate moves + self.test.run_test(toolhead, axis, gcmd) + raw_values = [(chip_axis, chip.finish_measurements()) + for chip_axis, chip in self.accel_chips + if axis in chip_axis or chip_axis in axis] + for chip_axis, chip_values in raw_values: + gcmd.respond_info("%s-axis accelerometer stats: %s" % ( + chip_axis, chip_values.get_stats(),)) + if not chip_values: + raise gcmd.error( + "%s-axis accelerometer measured no data" % ( + chip_axis,)) + new_data = helper.process_accelerometer_data(chip_values) + if calibration_data[axis] is None: + calibration_data[axis] = new_data + else: + calibration_data[axis].join(new_data) + + configfile = self.printer.lookup_object('configfile') + + for axis in calibrate_axes: + gcmd.respond_info( + "Calculating the best input shaper parameters for %s axis" + % (axis,)) + calibration_data[axis].normalize_to_frequencies() + shaper_name, shaper_freq, shapers_vals = helper.find_best_shaper( + calibration_data[axis], gcmd.respond_info) + gcmd.respond_info( + "Recommended shaper_type_%s = %s, shaper_freq_%s = %.1f Hz" + % (axis, shaper_name, axis, shaper_freq)) + helper.save_params(configfile, axis, shaper_name, shaper_freq) + csv_name = self.save_calibration_data( + 'calibration_data', name_suffix, helper, axis, + calibration_data[axis], shapers_vals) + gcmd.respond_info( + "Shaper calibration data written to %s file" % (csv_name,)) + + gcmd.respond_info( + "The SAVE_CONFIG command will update the printer config file\n" + "with these parameters and restart the printer.") + if input_shaper is not None: + input_shaper.enable_shaping() + gcmd.respond_info("Re-enabled [input_shaper] after calibration") + + def cmd_MEASURE_AXES_NOISE(self, gcmd): + meas_time = gcmd.get_float("MEAS_TIME", 2.) + for _, chip in self.accel_chips: + chip.start_measurements() + self.printer.lookup_object('toolhead').dwell(meas_time) + raw_values = [(axis, chip.finish_measurements()) + for axis, chip in self.accel_chips] + helper = shaper_calibrate.ShaperCalibrate(self.printer) + for axis, raw_data in raw_values: + data = helper.process_accelerometer_data(raw_data) + vx = data.psd_x.mean() + vy = data.psd_y.mean() + vz = data.psd_z.mean() + gcmd.respond_info("Axes noise for %s-axis accelerometer: " + "%.6f (x), %.6f (y), %.6f (z)" % ( + axis, vx, vy, vz)) + + def is_valid_name_suffix(self, name_suffix): + return name_suffix.replace('-', '').replace('_', '').isalnum() + + def get_filename(self, base, name_suffix, axis=None, point=None): + name = base + if axis: + name += '_' + axis + if point: + name += "_%.3f_%.3f_%.3f" % (point[0], point[1], point[2]) + name += '_' + name_suffix + return os.path.join("/tmp", name + ".csv") + + def save_calibration_data(self, base_name, name_suffix, shaper_calibrate, + axis, calibration_data, shapers_vals=None): + output = self.get_filename(base_name, name_suffix, axis) + shaper_calibrate.save_calibration_data(output, calibration_data, + shapers_vals) + return output + +def load_config(config): + return ResonanceTester(config) diff --git a/klippy/extras/shaper_calibrate.py b/klippy/extras/shaper_calibrate.py new file mode 100644 index 00000000..67160b4f --- /dev/null +++ b/klippy/extras/shaper_calibrate.py @@ -0,0 +1,382 @@ +# Automatic calibration of input shapers +# +# Copyright (C) 2020 Dmitry Butyugin +# +# This file may be distributed under the terms of the GNU GPLv3 license. +import importlib, logging, math, multiprocessing + +MIN_FREQ = 5. +MAX_FREQ = 200. +WINDOW_T_SEC = 0.5 +MAX_SHAPER_FREQ = 150. + +TEST_DAMPING_RATIOS=[0.075, 0.1, 0.15] +SHAPER_DAMPING_RATIO = 0.1 + +###################################################################### +# Input shapers +###################################################################### + +class InputShaperCfg: + def __init__(self, name, init_func, min_freq): + self.name = name + self.init_func = init_func + self.min_freq = min_freq + +def get_zv_shaper(shaper_freq, damping_ratio): + df = math.sqrt(1. - damping_ratio**2) + K = math.exp(-damping_ratio * math.pi / df) + t_d = 1. / (shaper_freq * df) + A = [1., K] + T = [0., .5*t_d] + return (A, T) + +def get_zvd_shaper(shaper_freq, damping_ratio): + df = math.sqrt(1. - damping_ratio**2) + K = math.exp(-damping_ratio * math.pi / df) + t_d = 1. / (shaper_freq * df) + A = [1., 2.*K, K**2] + T = [0., .5*t_d, t_d] + return (A, T) + +def get_mzv_shaper(shaper_freq, damping_ratio): + df = math.sqrt(1. - damping_ratio**2) + K = math.exp(-.75 * damping_ratio * math.pi / df) + t_d = 1. / (shaper_freq * df) + + a1 = 1. - 1. / math.sqrt(2.) + a2 = (math.sqrt(2.) - 1.) * K + a3 = a1 * K * K + + A = [a1, a2, a3] + T = [0., .375*t_d, .75*t_d] + return (A, T) + +def get_ei_shaper(shaper_freq, damping_ratio): + v_tol = 0.05 # vibration tolerance + df = math.sqrt(1. - damping_ratio**2) + K = math.exp(-damping_ratio * math.pi / df) + t_d = 1. / (shaper_freq * df) + + a1 = .25 * (1. + v_tol) + a2 = .5 * (1. - v_tol) * K + a3 = a1 * K * K + + A = [a1, a2, a3] + T = [0., .5*t_d, t_d] + return (A, T) + +def get_2hump_ei_shaper(shaper_freq, damping_ratio): + v_tol = 0.05 # vibration tolerance + df = math.sqrt(1. - damping_ratio**2) + K = math.exp(-damping_ratio * math.pi / df) + t_d = 1. / (shaper_freq * df) + + V2 = v_tol**2 + X = pow(V2 * (math.sqrt(1. - V2) + 1.), 1./3.) + a1 = (3.*X*X + 2.*X + 3.*V2) / (16.*X) + a2 = (.5 - a1) * K + a3 = a2 * K + a4 = a1 * K * K * K + + A = [a1, a2, a3, a4] + T = [0., .5*t_d, t_d, 1.5*t_d] + return (A, T) + +def get_3hump_ei_shaper(shaper_freq, damping_ratio): + v_tol = 0.05 # vibration tolerance + df = math.sqrt(1. - damping_ratio**2) + K = math.exp(-damping_ratio * math.pi / df) + t_d = 1. / (shaper_freq * df) + + K2 = K*K + a1 = 0.0625 * (1. + 3. * v_tol + 2. * math.sqrt(2. * (v_tol + 1.) * v_tol)) + a2 = 0.25 * (1. - v_tol) * K + a3 = (0.5 * (1. + v_tol) - 2. * a1) * K2 + a4 = a2 * K2 + a5 = a1 * K2 * K2 + + A = [a1, a2, a3, a4, a5] + T = [0., .5*t_d, t_d, 1.5*t_d, 2.*t_d] + return (A, T) + +INPUT_SHAPERS = [ + InputShaperCfg('zv', get_zv_shaper, 15.), + InputShaperCfg('mzv', get_mzv_shaper, 25.), + InputShaperCfg('ei', get_ei_shaper, 30.), + InputShaperCfg('2hump_ei', get_2hump_ei_shaper, 37.5), + InputShaperCfg('3hump_ei', get_3hump_ei_shaper, 50.), +] + +###################################################################### +# Frequency response calculation and shaper auto-tuning +###################################################################### + +class CalibrationData: + def __init__(self, freq_bins, psd_sum, psd_x, psd_y, psd_z): + self.freq_bins = freq_bins + self.psd_sum = psd_sum + self.psd_x = psd_x + self.psd_y = psd_y + self.psd_z = psd_z + self._psd_list = [self.psd_sum, self.psd_x, self.psd_y, self.psd_z] + self.data_sets = 1 + def join(self, other): + np = self.numpy + joined_data_sets = self.data_sets + other.data_sets + for psd, other_psd in zip(self._psd_list, other._psd_list): + # `other` data may be defined at different frequency bins, + # interpolating to fix that. + other_normalized = other.data_sets * np.interp( + self.freq_bins, other.freq_bins, other_psd) + psd *= self.data_sets + psd[:] = (psd + other_normalized) * (1. / joined_data_sets) + self.data_sets = joined_data_sets + def set_numpy(self, numpy): + self.numpy = numpy + def normalize_to_frequencies(self): + for psd in self._psd_list: + # Avoid division by zero errors + psd /= self.freq_bins + .1 + # Remove low-frequency noise + psd[self.freq_bins < MIN_FREQ] = 0. + + +class ShaperCalibrate: + def __init__(self, printer): + self.printer = printer + self.error = printer.command_error if printer else Exception + try: + self.numpy = importlib.import_module('numpy') + except ImportError: + raise self.error( + "Failed to import `numpy` module, make sure it was " + "installed via `~/klippy-env/bin/pip install` (refer to " + "docs/Measuring_Resonances.md for more details).") + + def background_process_exec(self, method, args): + if self.printer is None: + return method(*args) + import queuelogger + parent_conn, child_conn = multiprocessing.Pipe() + def wrapper(): + queuelogger.clear_bg_logging() + try: + res = method(*args) + except: + child_conn.send((True, traceback.format_exc())) + child_conn.close() + return + child_conn.send((False, res)) + child_conn.close() + # Start a process to perform the calculation + calc_proc = multiprocessing.Process(target=wrapper) + calc_proc.daemon = True + calc_proc.start() + # Wait for the process to finish + reactor = self.printer.get_reactor() + gcode = self.printer.lookup_object("gcode") + eventtime = last_report_time = reactor.monotonic() + while calc_proc.is_alive(): + if eventtime > last_report_time + 5.: + last_report_time = eventtime + gcode.respond_info("Wait for calculations..", log=False) + eventtime = reactor.pause(eventtime + .1) + # Return results + is_err, res = parent_conn.recv() + if is_err: + raise self.error("Error in remote calculation: %s" % (res,)) + calc_proc.join() + parent_conn.close() + return res + + def _split_into_windows(self, x, window_size, overlap): + # Memory-efficient algorithm to split an input 'x' into a series + # of overlapping windows + step_between_windows = window_size - overlap + n_windows = (x.shape[-1] - overlap) // step_between_windows + shape = (window_size, n_windows) + strides = (x.strides[-1], step_between_windows * x.strides[-1]) + return self.numpy.lib.stride_tricks.as_strided( + x, shape=shape, strides=strides, writeable=False) + + def _psd(self, x, fs, nfft): + # Calculate power spectral density (PSD) using Welch's algorithm + np = self.numpy + window = np.kaiser(nfft, 6.) + # Compensation for windowing loss + scale = 1.0 / (window**2).sum() + + # Split into overlapping windows of size nfft + overlap = nfft // 2 + x = self._split_into_windows(x, nfft, overlap) + + # First detrend, then apply windowing function + x = window[:, None] * (x - np.mean(x, axis=0)) + + # Calculate frequency response for each window using FFT + result = np.fft.rfft(x, n=nfft, axis=0) + result = np.conjugate(result) * result + result *= scale / fs + # For one-sided FFT output the response must be doubled, except + # the last point for unpaired Nyquist frequency (assuming even nfft) + # and the 'DC' term (0 Hz) + result[1:-1,:] *= 2. + + # Welch's algorithm: average response over windows + psd = result.real.mean(axis=-1) + + # Calculate the frequency bins + freqs = np.fft.rfftfreq(nfft, 1. / fs) + return freqs, psd + + def calc_freq_response(self, raw_values): + np = self.numpy + if raw_values is None: + return None + if isinstance(raw_values, np.ndarray): + data = raw_values + else: + data = np.array(raw_values.decode_samples()) + + N = data.shape[0] + T = data[-1,0] - data[0,0] + SAMPLING_FREQ = N / T + # Round up to the nearest power of 2 for faster FFT + M = 1 << int(SAMPLING_FREQ * WINDOW_T_SEC - 1).bit_length() + if N <= M: + return None + + # Calculate PSD (power spectral density) of vibrations per + # frequency bins (the same bins for X, Y, and Z) + fx, px = self._psd(data[:,1], SAMPLING_FREQ, M) + fy, py = self._psd(data[:,2], SAMPLING_FREQ, M) + fz, pz = self._psd(data[:,3], SAMPLING_FREQ, M) + return CalibrationData(fx, px+py+pz, px, py, pz) + + def process_accelerometer_data(self, data): + calibration_data = self.background_process_exec( + self.calc_freq_response, (data,)) + if calibration_data is None: + raise self.error( + "Internal error processing accelerometer data %s" % (data,)) + calibration_data.set_numpy(self.numpy) + return calibration_data + + def _estimate_shaper(self, shaper, test_damping_ratio, test_freqs): + np = self.numpy + + A, T = np.array(shaper[0]), np.array(shaper[1]) + inv_D = 1. / A.sum() + + omega = 2. * math.pi * test_freqs + damping = test_damping_ratio * omega + omega_d = omega * math.sqrt(1. - test_damping_ratio**2) + W = A * np.exp(np.outer(-damping, (T[-1] - T))) + S = W * np.sin(np.outer(omega_d, T)) + C = W * np.cos(np.outer(omega_d, T)) + return np.sqrt(S.sum(axis=1)**2 + C.sum(axis=1)**2) * inv_D + + def _estimate_remaining_vibrations(self, shaper, test_damping_ratio, + freq_bins, psd): + vals = self._estimate_shaper(shaper, test_damping_ratio, freq_bins) + remaining_vibrations = (vals * psd).sum() / psd.sum() + return (remaining_vibrations, vals) + + def fit_shaper(self, shaper_cfg, calibration_data): + np = self.numpy + + test_freqs = np.arange(shaper_cfg.min_freq, MAX_SHAPER_FREQ, .2) + + freq_bins = calibration_data.freq_bins + psd = calibration_data.psd_sum[freq_bins <= MAX_FREQ] + freq_bins = freq_bins[freq_bins <= MAX_FREQ] + + best_freq = None + best_remaining_vibrations = 0 + best_shaper_vals = [] + + for test_freq in test_freqs[::-1]: + cur_remaining_vibrations = 0. + shaper_vals = np.zeros(shape=freq_bins.shape) + shaper = shaper_cfg.init_func(test_freq, SHAPER_DAMPING_RATIO) + # Exact damping ratio of the printer is unknown, pessimizing + # remaining vibrations over possible damping values. + for dr in TEST_DAMPING_RATIOS: + vibrations, vals = self._estimate_remaining_vibrations( + shaper, dr, freq_bins, psd) + shaper_vals = np.maximum(shaper_vals, vals) + if vibrations > cur_remaining_vibrations: + cur_remaining_vibrations = vibrations + if (best_freq is None or + best_remaining_vibrations > cur_remaining_vibrations): + # The current frequency is better for the shaper. + best_freq = test_freq + best_remaining_vibrations = cur_remaining_vibrations + best_shaper_vals = shaper_vals + return (best_freq, best_remaining_vibrations, best_shaper_vals) + + def find_best_shaper(self, calibration_data, logger=None): + best_shaper = prev_shaper = None + best_freq = prev_freq = 0. + best_vibrations = prev_vibrations = 0. + all_shaper_vals = [] + for shaper in INPUT_SHAPERS: + shaper_freq, vibrations, shaper_vals = self.background_process_exec( + self.fit_shaper, (shaper, calibration_data)) + if logger is not None: + logger("Fitted shaper '%s' frequency = %.1f Hz " + "(vibrations = %.1f%%)" % ( + shaper.name, shaper_freq, vibrations * 100.)) + if best_shaper is None or 1.75 * vibrations < best_vibrations: + if 1.25 * vibrations < prev_vibrations: + best_shaper = shaper.name + best_freq = shaper_freq + best_vibrations = vibrations + else: + # The current shaper is good, but not sufficiently better + # than the previous one, using previous shaper instead. + best_shaper = prev_shaper + best_freq = prev_freq + best_vibrations = prev_vibrations + prev_shaper = shaper.name + prev_shaper_vals = shaper_vals + prev_freq = shaper_freq + prev_vibrations = vibrations + all_shaper_vals.append((shaper.name, shaper_freq, shaper_vals)) + return (best_shaper, best_freq, all_shaper_vals) + + def save_params(self, configfile, axis, shaper_name, shaper_freq): + if axis == 'xy': + self.save_params(configfile, 'x', shaper_name, shaper_freq) + self.save_params(configfile, 'y', shaper_name, shaper_freq) + else: + configfile.set('input_shaper', 'shaper_type_'+axis, shaper_name) + configfile.set('input_shaper', 'shaper_freq_'+axis, + '%.1f' % (shaper_freq,)) + + def save_calibration_data(self, output, calibration_data, + shapers_vals=None): + try: + with open(output, "w") as csvfile: + csvfile.write("freq,psd_x,psd_y,psd_z,psd_xyz") + if shapers_vals: + for name, freq, _ in shapers_vals: + csvfile.write(",%s(%.1f)" % (name, freq)) + csvfile.write("\n") + num_freqs = calibration_data.freq_bins.shape[0] + for i in range(num_freqs): + if calibration_data.freq_bins[i] >= MAX_FREQ: + break + csvfile.write("%.1f,%.3e,%.3e,%.3e,%.3e" % ( + calibration_data.freq_bins[i], + calibration_data.psd_x[i], + calibration_data.psd_y[i], + calibration_data.psd_z[i], + calibration_data.psd_sum[i])) + if shapers_vals: + for _, _, vals in shapers_vals: + csvfile.write(",%.3f" % (vals[i],)) + csvfile.write("\n") + except IOError as e: + raise self.error("Error writing to file '%s': %s", output, str(e)) diff --git a/scripts/calibrate_shaper.py b/scripts/calibrate_shaper.py new file mode 100755 index 00000000..22b32b23 --- /dev/null +++ b/scripts/calibrate_shaper.py @@ -0,0 +1,166 @@ +#!/usr/bin/env python2 +# Shaper auto-calibration script +# +# Copyright (C) 2020 Dmitry Butyugin +# Copyright (C) 2020 Kevin O'Connor +# +# This file may be distributed under the terms of the GNU GPLv3 license. +from __future__ import print_function +import optparse, os, sys +import numpy as np, matplotlib +sys.path.append(os.path.join(os.path.dirname(os.path.realpath(__file__)), + '..', 'klippy', 'extras')) +from shaper_calibrate import CalibrationData, ShaperCalibrate + +def parse_log(logname): + with open(logname) as f: + for header in f: + if not header.startswith('#'): + break + if not header.startswith('freq,psd_x,psd_y,psd_z,psd_xyz'): + # Raw accelerometer data + return np.loadtxt(logname, comments='#', delimiter=',') + # Parse power spectral density data + data = np.loadtxt(logname, skiprows=1, comments='#', delimiter=',') + calibration_data = CalibrationData( + freq_bins=data[:,0], psd_sum=data[:,4], + psd_x=data[:,1], psd_y=data[:,2], psd_z=data[:,3]) + calibration_data.set_numpy(np) + # If input shapers are present in the CSV file, the frequency + # response is already normalized to input frequencies + if 'mzv' not in header: + calibration_data.normalize_to_frequencies() + return calibration_data + +###################################################################### +# Shaper calibration +###################################################################### + +# Find the best shaper parameters +def calibrate_shaper(datas, csv_output): + helper = ShaperCalibrate(printer=None) + if isinstance(datas[0], CalibrationData): + calibration_data = datas[0] + for data in datas[1:]: + calibration_data.join(data) + else: + # Process accelerometer data + calibration_data = helper.process_accelerometer_data(datas[0]) + for data in datas[1:]: + calibration_data.join(helper.process_accelerometer_data(data)) + calibration_data.normalize_to_frequencies() + shaper_name, shaper_freq, shapers_vals = helper.find_best_shaper( + calibration_data, print) + print("Recommended shaper is %s @ %.1f Hz" % (shaper_name, shaper_freq)) + if csv_output is not None: + helper.save_calibration_data( + csv_output, calibration_data, shapers_vals) + return shaper_name, shapers_vals, calibration_data + +###################################################################### +# Plot frequency response and suggested input shapers +###################################################################### + +def plot_freq_response(calibration_data, shapers_vals, + selected_shaper, max_freq): + freqs = calibration_data.freq_bins + psd = calibration_data.psd_sum[freqs <= max_freq] + px = calibration_data.psd_x[freqs <= max_freq] + py = calibration_data.psd_y[freqs <= max_freq] + pz = calibration_data.psd_z[freqs <= max_freq] + freqs = freqs[freqs <= max_freq] + + fontP = matplotlib.font_manager.FontProperties() + fontP.set_size('x-small') + + fig, ax = matplotlib.pyplot.subplots() + ax.set_xlabel('Frequency, Hz') + ax.set_xlim([0, max_freq]) + ax.set_ylabel('Power spectral density') + + ax.plot(freqs, psd, label='X+Y+Z', color='purple') + ax.plot(freqs, px, label='X', color='red') + ax.plot(freqs, py, label='Y', color='green') + ax.plot(freqs, pz, label='Z', color='blue') + + if shapers_vals: + ax.set_title("Frequency response and shapers") + else: + ax.set_title("Frequency response") + ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator()) + ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator()) + ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator()) + ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator()) + ax.ticklabel_format(axis='y', style='scientific', scilimits=(0,0)) + ax.grid(which='major', color='grey') + ax.grid(which='minor', color='lightgrey') + ax.legend(loc='upper right', prop=fontP) + + if shapers_vals: + ax2 = ax.twinx() + ax2.set_ylabel('Shaper vibration reduction (ratio)') + best_shaper_vals = None + for name, freq, vals in shapers_vals: + label = "%s (%.1f Hz)" % (name.upper(), freq) + linestyle = 'dotted' + if name == selected_shaper: + label += ' (selected)' + linestyle = 'dashdot' + best_shaper_vals = vals + ax2.plot(freqs, vals, label=label, linestyle=linestyle) + ax.plot(freqs, psd * best_shaper_vals, + label='After\nshaper', color='cyan') + ax2.legend(loc='upper left', prop=fontP) + + fig.tight_layout() + return fig + +###################################################################### +# Startup +###################################################################### + +def setup_matplotlib(output_to_file): + global matplotlib + if output_to_file: + matplotlib.rcParams.update({'figure.autolayout': True}) + matplotlib.use('Agg') + import matplotlib.pyplot, matplotlib.dates, matplotlib.font_manager + import matplotlib.ticker + +def main(): + # Parse command-line arguments + usage = "%prog [options] " + opts = optparse.OptionParser(usage) + opts.add_option("-o", "--output", type="string", dest="output", + default=None, help="filename of output graph") + opts.add_option("-c", "--csv", type="string", dest="csv", + default=None, help="filename of output csv file") + opts.add_option("-f", "--max_freq", type="float", default=200., + help="maximum frequency to graph") + options, args = opts.parse_args() + if len(args) < 1: + opts.error("Incorrect number of arguments") + + # Parse data + datas = [parse_log(fn) for fn in args] + + # Calibrate shaper and generate outputs + selected_shaper, shapers_vals, calibration_data = calibrate_shaper( + datas, options.csv) + + if not options.csv or options.output: + # Draw graph + setup_matplotlib(options.output is not None) + + fig = plot_freq_response(calibration_data, shapers_vals, + selected_shaper, options.max_freq) + + # Show graph + if options.output is None: + matplotlib.pyplot.show() + else: + fig.set_size_inches(8, 6) + fig.savefig(options.output) + +if __name__ == '__main__': + main() diff --git a/scripts/graph_accelerometer.py b/scripts/graph_accelerometer.py index e544d0f7..0e975bf5 100755 --- a/scripts/graph_accelerometer.py +++ b/scripts/graph_accelerometer.py @@ -2,38 +2,35 @@ # Generate adxl345 accelerometer graphs # # Copyright (C) 2020 Kevin O'Connor +# Copyright (C) 2020 Dmitry Butyugin # # This file may be distributed under the terms of the GNU GPLv3 license. -import optparse -import matplotlib +import optparse, os, sys +from textwrap import wrap +import numpy as np, matplotlib +sys.path.append(os.path.join(os.path.dirname(os.path.realpath(__file__)), + '..', 'klippy', 'extras')) +from shaper_calibrate import ShaperCalibrate + +MAX_TITLE_LENGTH=80 def parse_log(logname): - f = open(logname, 'r') - out = [] - for line in f: - if line.startswith('#'): - continue - parts = line.split(',') - if len(parts) != 4: - continue - try: - fparts = [float(p) for p in parts] - except ValueError: - continue - out.append(fparts) - return out + return np.loadtxt(logname, comments='#', delimiter=',') + +###################################################################### +# Raw accelerometer graphing +###################################################################### def plot_accel(data, logname): - half_smooth_samples = 15 - first_time = data[0][0] - times = [d[0] - first_time for d in data] + first_time = data[0, 0] + times = data[:,0] - first_time fig, axes = matplotlib.pyplot.subplots(nrows=3, sharex=True) - axes[0].set_title("Accelerometer data (%s)" % (logname,)) + axes[0].set_title("\n".join(wrap("Accelerometer data (%s)" % (logname,), + MAX_TITLE_LENGTH))) axis_names = ['x', 'y', 'z'] for i in range(len(axis_names)): - adata = [d[i+1] for d in data] - avg = sum(adata) / len(adata) - adata = [ad - avg for ad in adata] + avg = data[:,i+1].mean() + adata = data[:,i+1] - data[:,i+1].mean() ax = axes[i] ax.plot(times, adata, alpha=0.8) ax.grid(True) @@ -42,9 +39,145 @@ def plot_accel(data, logname): fig.tight_layout() return fig -def setup_matplotlib(output_to_file): + +###################################################################### +# Frequency graphing +###################################################################### + +# Calculate estimated "power spectral density" +def calc_freq_response(data, max_freq): + helper = ShaperCalibrate(printer=None) + return helper.process_accelerometer_data(data) + +def calc_specgram(data, axis): + N = data.shape[0] + Fs = N / (data[-1,0] - data[0,0]) + # Round up to a power of 2 for faster FFT + M = 1 << int(.5 * Fs - 1).bit_length() + window = np.kaiser(M, 6.) + def _specgram(x): + return matplotlib.mlab.specgram( + x, Fs=Fs, NFFT=M, noverlap=M//2, window=window, + mode='psd', detrend='mean', scale_by_freq=False) + + d = {'x': data[:,1], 'y': data[:,2], 'z': data[:,3]} + if axis != 'all': + pdata, bins, t = _specgram(d[axis]) + else: + pdata, bins, t = _specgram(d['x']) + for ax in 'yz': + pdata += _specgram(d[ax])[0] + return pdata, bins, t + +def plot_frequency(datas, lognames, max_freq): + calibration_data = calc_freq_response(datas[0], max_freq) + for data in datas[1:]: + calibration_data.join(calc_freq_response(data, max_freq)) + freqs = calibration_data.freq_bins + psd = calibration_data.psd_sum[freqs <= max_freq] + px = calibration_data.psd_x[freqs <= max_freq] + py = calibration_data.psd_y[freqs <= max_freq] + pz = calibration_data.psd_z[freqs <= max_freq] + freqs = freqs[freqs <= max_freq] + + fig, ax = matplotlib.pyplot.subplots() + ax.set_title("\n".join(wrap( + "Frequency response (%s)" % (', '.join(lognames)), MAX_TITLE_LENGTH))) + ax.set_xlabel('Frequency (Hz)') + ax.set_ylabel('Power spectral density') + + ax.plot(freqs, psd, label='X+Y+Z', alpha=0.6) + ax.plot(freqs, px, label='X', alpha=0.6) + ax.plot(freqs, py, label='Y', alpha=0.6) + ax.plot(freqs, pz, label='Z', alpha=0.6) + + ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator()) + ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator()) + ax.grid(which='major', color='grey') + ax.grid(which='minor', color='lightgrey') + ax.ticklabel_format(axis='y', style='scientific', scilimits=(0,0)) + + fontP = matplotlib.font_manager.FontProperties() + fontP.set_size('x-small') + ax.legend(loc='best', prop=fontP) + fig.tight_layout() + return fig + +def plot_compare_frequency(datas, lognames, max_freq): + fig, ax = matplotlib.pyplot.subplots() + ax.set_title('Frequency responses comparison') + ax.set_xlabel('Frequency (Hz)') + ax.set_ylabel('Power spectral density') + + for data, logname in zip(datas, lognames): + calibration_data = calc_freq_response(data, max_freq) + freqs = calibration_data.freq_bins + psd = calibration_data.psd_sum[freqs <= max_freq] + freqs = freqs[freqs <= max_freq] + ax.plot(freqs, psd, label="\n".join(wrap(logname, 60)), alpha=0.6) + + ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator()) + ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator()) + ax.grid(which='major', color='grey') + ax.grid(which='minor', color='lightgrey') + fontP = matplotlib.font_manager.FontProperties() + fontP.set_size('x-small') + ax.legend(loc='best', prop=fontP) + fig.tight_layout() + return fig + +# Plot data in a "spectrogram colormap" +def plot_specgram(data, logname, max_freq, axis): + pdata, bins, t = calc_specgram(data, axis) + + fig, ax = matplotlib.pyplot.subplots() + ax.set_title("\n".join(wrap("Spectrogram %s (%s)" % (axis, logname), + MAX_TITLE_LENGTH))) + ax.pcolormesh(t, bins, pdata, norm=matplotlib.colors.LogNorm()) + ax.set_ylim([0., max_freq]) + ax.set_ylabel('frequency (hz)') + ax.set_xlabel('Time (s)') + fig.tight_layout() + return fig + +###################################################################### +# CSV output +###################################################################### + +def write_frequency_response(datas, output): + helper = ShaperCalibrate(printer=None) + calibration_data = helper.process_accelerometer_data(datas[0]) + for data in datas[1:]: + calibration_data.join(helper.process_accelerometer_data(data)) + helper.save_calibration_data(output, calibration_data) + +def write_specgram(psd, freq_bins, time, output): + M = freq_bins.shape[0] + with open(output, "w") as csvfile: + csvfile.write("freq\\t") + for ts in time: + csvfile.write(",%.6f" % (ts,)) + csvfile.write("\n") + for i in range(M): + csvfile.write("%.1f" % (freq_bins[i],)) + for value in psd[i,:]: + csvfile.write(",%.6e" % (value,)) + csvfile.write("\n") + +###################################################################### +# Startup +###################################################################### + +def is_csv_output(output): + return output and os.path.splitext(output)[1].lower() == '.csv' + +def setup_matplotlib(output): global matplotlib - if output_to_file: + if is_csv_output(output): + # Only mlab may be necessary with CSV output + import matplotlib.mlab + return + if output: matplotlib.rcParams.update({'figure.autolayout': True}) matplotlib.use('Agg') import matplotlib.pyplot, matplotlib.dates, matplotlib.font_manager @@ -52,20 +185,51 @@ def setup_matplotlib(output_to_file): def main(): # Parse command-line arguments - usage = "%prog [options] " + usage = "%prog [options] " opts = optparse.OptionParser(usage) opts.add_option("-o", "--output", type="string", dest="output", default=None, help="filename of output graph") + opts.add_option("-f", "--max_freq", type="float", default=200., + help="maximum frequency to graph") + opts.add_option("-r", "--raw", action="store_true", + help="graph raw accelerometer data") + opts.add_option("-c", "--compare", action="store_true", + help="graph comparison of power spectral density " + "between different accelerometer data files") + opts.add_option("-s", "--specgram", action="store_true", + help="graph spectrogram of accelerometer data") + opts.add_option("-a", type="string", dest="axis", default="all", + help="axis to graph (one of 'all', 'x', 'y', or 'z')") options, args = opts.parse_args() - if len(args) != 1: + if len(args) < 1: opts.error("Incorrect number of arguments") # Parse data - data = parse_log(args[0]) + datas = [parse_log(fn) for fn in args] + + setup_matplotlib(options.output) + + if is_csv_output(options.output): + if options.raw: + opts.error("raw mode is not supported with csv output") + if options.compare: + opts.error("comparison mode is not supported with csv output") + if options.specgram: + pdata, bins, t = calc_specgram(datas[0], options.axis) + write_specgram(pdata, bins, t, options.output) + else: + write_frequency_response(datas, options.output) + return # Draw graph - setup_matplotlib(options.output is not None) - fig = plot_accel(data, args[0]) + if options.raw: + fig = plot_accel(datas[0], args[0]) + elif options.specgram: + fig = plot_specgram(datas[0], args[0], options.max_freq, options.axis) + elif options.compare: + fig = plot_compare_frequency(datas, args, options.max_freq) + else: + fig = plot_frequency(datas, args, options.max_freq) # Show graph if options.output is None: