10 KiB
This notebook analyzes a video capture from the laser camera where a far mirror was set up 2.35m across the room. The laser diode and camera were both placed on the vibration isolated mount of the STM. In front of the laser diode, a ø2mm apterture was placed to somewhat contain the beam spread.
From the video, we recover the center position of the laser in each frame. Then we can analyze the movement of the laser beam over time and extract frequency content. This should give us an indication of the transmittance of the vibration isolation over time.
import cv2
import numpy as np
import matplotlib.pyplot as plt
from scipy import ndimage
import time
capture = cv2.VideoCapture("./laser-img/2025-09-07-011400.webm")
Fetch video capture metadata for processing:
frame_count = int(capture.get(cv2.CAP_PROP_FRAME_COUNT))
frame_rate = int(capture.get(cv2.CAP_PROP_FPS))
frame_time = 1 / frame_rate
frame_width = int(capture.get(cv2.CAP_PROP_FRAME_WIDTH))
frame_height = int(capture.get(cv2.CAP_PROP_FRAME_HEIGHT))
# Overwrite frame_count for limiting analysis window
frame_count = 3000
print(f"Capture Info: {frame_count} frames with dimension {frame_width}×{frame_height} at {frame_rate} Hz, totalling {frame_count * frame_time:.2f} seconds")
Capture Info: 3000 frames with dimension 1920×1080 at 30 Hz, totalling 100.00 seconds
Laser Beam Position Recovery
The laser beam position is recovered using the ndimage.center_of_mass() function. This function hinges on the background of the image being completely dark — otherwise there will be a positional bias towards the center of the image.
def find_laser_beam(image):
image_monochrome = np.mean(image, axis=2)
intensity = np.mean(image)
y, x = ndimage.center_of_mass(image_monochrome)
return np.array((x, y, intensity))
Video Capture Processing
Process the capture frame by frame to recover the laser beam position in each frame. For verification of the recovery function, a few shots from the capture are visualized.
processing_start = time.monotonic()
laser_position = np.empty((frame_count, 3), np.dtype('float32'))
current_frame = 0
ret = True
while current_frame < frame_count and ret:
ret, frame = capture.read()
if ret:
current_laser = find_laser_beam(frame)
current_frame_time = current_frame * frame_time
laser_position[current_frame] = current_laser
if current_frame % (frame_count // 9) == 0:
print(f"Snapshot capture of frame {current_frame:6} at {current_frame_time:6.2f}s: Intensity {current_laser[2]:.3f}")
plt.figure()
plt.imshow(frame)
plt.axvline(x=current_laser[0], color="red")
plt.axhline(y=current_laser[1], color="red")
plt.figtext(0.5, 0.05, f"Frame {current_frame} at {current_frame_time:.2f} seconds; Intensity {current_laser[2]:.3f}", ha="center")
current_frame += 1
processing_time = time.monotonic() - processing_start
print(f"Video processing took {processing_time:.2f} seconds ({processing_time / frame_count:.3f} seconds per frame).")
Snapshot capture of frame 0 at 0.00s: Intensity 9.282
Snapshot capture of frame 333 at 11.10s: Intensity 10.440
Snapshot capture of frame 666 at 22.20s: Intensity 12.150
Snapshot capture of frame 999 at 33.30s: Intensity 12.801
Snapshot capture of frame 1332 at 44.40s: Intensity 12.652
Snapshot capture of frame 1665 at 55.50s: Intensity 11.099
Snapshot capture of frame 1998 at 66.60s: Intensity 11.897
Snapshot capture of frame 2331 at 77.70s: Intensity 12.612
Snapshot capture of frame 2664 at 88.80s: Intensity 12.257
Snapshot capture of frame 2997 at 99.90s: Intensity 12.728
Video processing took 300.84 seconds (0.100 seconds per frame).
Derive data for the plots below from the raw analysis.
time_axis = np.linspace(0, frame_count * frame_time, frame_count)
laser_displacement_zeroed = laser_position[:, 0:2] - np.nanmean(laser_position[:, 0:2], axis=0)
velocity_xy = np.diff(laser_position[:, 0:2])
velocity_abs = np.nan_to_num(np.linalg.norm(velocity_xy, axis=1))
max_intensity = np.max(laser_position[:, 2])
laser_intensity = laser_position[:, 2] / max_intensity
Time-dependent Laser Position and Intensity plots
The following figure shows the laser displacement, velocity, and intensity over the whole capture.
fig, ax = plt.subplots(4, 1, figsize=(18, 16))
ax_displacement = ax[0:3]
ax_intensity = ax[3]
ax_displacement[0].plot(time_axis, laser_displacement_zeroed[:, 0] / frame_width, label="Y position")
ax_displacement[1].plot(time_axis, laser_displacement_zeroed[:, 1] / frame_width, label="X position")
ax_displacement[2].plot(time_axis, velocity_abs, label="Velocity")
for ax_d in ax_displacement:
ax_d.set_xlabel("time / seconds")
ax_d.set_ylabel("displacment / sensor width")
ax_d.legend()
ax_intensity.plot(time_axis, laser_intensity, label="Intensity")
ax_intensity.set_xlabel("time / seconds")
ax_intensity.set_ylabel("relative intensity")
ax_intensity.set_ylim(0, 1)
ax_intensity.legend()
fig.text(0.5, 0.05, "Laser X/Y, Velocity, and Intensity over time", ha="center")
pass
from scipy import signal
frequency_x, power_spectral_density_x = signal.periodogram(np.nan_to_num(laser_position[:, 1]), frame_rate)
frequency_y, power_spectral_density_y = signal.periodogram(np.nan_to_num(laser_position[:, 0]), frame_rate)
amplitude_spectral_density_x = np.sqrt(power_spectral_density_x)
amplitude_spectral_density_y = np.sqrt(power_spectral_density_y)
Spectral Analysis
The following figure looks at the spectral density in the displacement along X and Y axis. This gives us a good insight to see which vibration frequencies were transmitted through the isolation.
fig, ax = plt.subplots(2, 1, figsize=(16, 8))
ax[1].axvline(x=2.822, color="red")
ax[1].annotate(" theoretical spring resonance", xy=(2.822, 1080 / 2 - 150))
ax[0].loglog(frequency_x, amplitude_spectral_density_x, label="X axis")
ax[1].loglog(frequency_y, amplitude_spectral_density_y, label="Y axis")
for a in ax:
a.set_xlabel("frequency / Hz")
a.set_ylabel("amplitude spectral density")
a.set_xlim(0.5, 15)
a.set_ylim(1, 1080 / 2)
a.legend()
a.grid(axis="x", which="both")
a.grid(axis="y", which="major")
a.axvline(x=1.317, color="red")
a.annotate(" pendulum mode", xy=(1.317, 1080 / 2 - 150))
pass
Of note is the highlighted theoretical spring resonance which we had previously calculated using the known spring constant and mass of the system. Seeing a very visible peak at this frequency along the Y axis (the axis of the springs) is very nice.
In addition, the pendulum swinging resonance of the assembly was also roughly calculated and drawn into the diagrams. A peak is visible near this frequency as well.
# If you wish to save the extracted displacement data, use this:
# np.save("./laser-results/2025-09-07-011400-analyzed.npy", laser_position)
Theoretical Transmissibility vs. Measurement
As a further step, the code below calculates the theoretical transmissibility curve and then plots it overlayed with the laser displacement amplitude spectral density in Y direction.
This gives a nice comparison, but the results should be interpreted with care: We do not have information on the spectrum of the excitation that was present during the measurement. For sure it was not even density broadband noise. Additionally, towards lower amplitudes, the measurement data gets increasingly noisy because the signal reaches the limit of the camera resolution.
import math
import numpy as np
from pint import UnitRegistry
unit = UnitRegistry()
unit.formatter.default_format = "~"
# Parameters
spring_constant = 1.1 * 4 * unit.N / unit.mm
spring_length_resting = 112 * unit.mm
weight_total = 14 * unit.kg
dampening = 1 * unit.N / (unit.m / unit.s)
def spring_length_at(weight):
return (weight * unit.standard_gravity / spring_constant + spring_length_resting).to(unit.mm)
def resonant_freq_at(weight):
return (1 / (2 * math.pi) * np.sqrt(spring_constant / weight)).to(unit.Hz)
spring_length = spring_length_at(weight_total)
f0 = resonant_freq_at(weight_total)
print(f"Length: {spring_length:~.1f}")
print(f"Freq: {f0:~.3f}")
def lehr_dampening_factor(d, k, m):
return d / (2 * np.sqrt(m * k))
lehr_dampening = lehr_dampening_factor(dampening, spring_constant, weight_total)
lehr_dampening.ito_reduced_units()
def amplitude_ratio(lehr, f0, f):
eta = f / f0
return 1 / np.sqrt((1 - eta**2)**2 + (2 * eta * lehr)**2)
f_in = np.geomspace(0.5, 90, 100) * unit.Hz
ratio = amplitude_ratio(lehr_dampening, f0, f_in)
Length: 143.2 mm
Freq: 2.822 Hz
fig, ax = plt.subplots()
ax.axvline(x=2.822, color="red")
# ax.annotate(" theoretical spring resonance", xy=(2.822, 1080 / 2 - 150))
ax.loglog(frequency_y, amplitude_spectral_density_y, label="laser")
ax.loglog(f_in, ratio * 8, label="transmissibility")
a = ax
a.set_xlabel("frequency / Hz")
# a.set_ylabel("amplitude spectral density")
a.set_xlim(0.5, 15)
a.set_ylim(0.1, 1080 / 2)
a.legend()
a.grid(axis="x", which="both")
a.grid(axis="y", which="major")
pass
One more note about this graph: As the input amplitude is unknown, the theoretical curve was roughly scaled to sit ontop of the measurement data (ratio * 8 in the code above). This is not a proper curve fit.
Also keep in mind that the measurement data gets pretty meaningless below an amplitude of 1 (pixel) due to the resolution limit of the camera.