校表仪 tg-timer 算法笔记

**节拍误差：**

![](/media/202412/image-20240826102408879_1735263067.png)

![](/media/202412/image-20240905145819162_1735263072.png)

![](/media/202412/image-20240905145901252_1735263082.png)

![摆轮升角](https://productimg.xbiao.com/baike/2012/1228/20121228135666554046324.jpg)

### Period周期计算

```c
// 峰值检测器
static int peak_detector(float *buff, int a, int b)
{
	int i_max;
	double max = vmax(buff, a, b+1, &i_max);
	if(max <= 0) return -1;

int i;
	float v[b-a+1];
	memcpy(v, buff + a, sizeof(v));
	quickselect(v, b-a+1, (b-a+1)/2);
	float med = v[(b-a+1)/2];

for(i=a+1; i<i_max; i++)
		if(buff[i] <= med) break;
	if(i==i_max) return -1;
	for(i=i_max+1; i<=b; i++)
		if(buff[i] <= med) break;
	if(i==b+1) return -1;

int cnt = 0, down = 1;
	for(i=a+1; i<=b; i++) {
		if(buff[i] > (max + med) / 2) {
			cnt += down;
			down = 0;
		}
		if(buff[i] < med)
			down = 1;
	}
	debug("max = %f med = %f cnt = %d\n",max,med,cnt);
	if(cnt > 20)
		return -1;
	return i_max;
}
```

```c
static int compute_period(struct processing_buffers *b, int bph)
{
	double estimate;
	if(bph)
		estimate = peak_detector(b->samples_sc,
				7200 * b->sample_rate / bph - b->sample_rate / 50,
				7200 * b->sample_rate / bph + b->sample_rate / 50);
	else
		estimate = estimate_period(b);
	if(estimate == -1) {
		debug("failed to estimate period\n");
		return 1;
	}
	double delta = b->sample_rate * 0.02;
	double new_estimate = estimate;
	double sum = 0;
	double sq_sum = 0;
	int count = 0;
	int cycle = 1;
	for(;;) {
		int inf = floor(new_estimate * cycle - delta);
		int sup = ceil(new_estimate * cycle + delta);
		if(sup > b->sample_count * 2 / 3)
			break;
		new_estimate = peak_detector(b->samples_sc,inf,sup);
		if(new_estimate == -1) {
			debug("cycle = %d peak not found\n",cycle);
			return 1;
		}
		new_estimate /= cycle;
		if(new_estimate < estimate - delta || new_estimate > estimate + delta) {
			debug("cycle = %d new_estimate = %f invalid peak\n",cycle,new_estimate/b->sample_rate);
			return 1;
		} else
			debug("cycle = %d new_estimate = %f\n",cycle,new_estimate/b->sample_rate);
		if(inf > b->sample_count / 3) {
			sum += new_estimate;
			sq_sum += new_estimate * new_estimate;
			count++;
		}
		cycle++;
	}
	if(count > 0) estimate = sum / count;
	b->period = estimate;
	if(count > 1)
		b->sigma = sqrt((sq_sum - count * estimate * estimate)/ (count-1));
	else
		b->sigma = 0;	// No std. dev. estimate possible with just 1 sample
	return 0;
}

```

### Phase相位计算

```c
static void compute_phase(struct processing_buffers *p, double period)
{
	int i;
	double x = 0, y = 0;
	for(i = 0; i < period; i++) {
		int j;
		p->waveform[i] = 0;
		for(j=0;;j++) {
			int n = round(i + j * period);
			if(n >= p->sample_count) break;
			p->waveform[i] += p->samples[n];
		}
		p->waveform[i] /= j;
	}
	for(i=0; i<period; i++) {
		double a = i * 2 * M_PI / period;
		x += p->waveform[i] * cos(a);
		y += p->waveform[i] * sin(a);
	}
	p->phase = period * (M_PI + atan2(y,x)) / (2 * M_PI);
}
```

###  Wave波形计算

```c
static void compute_waveform(struct processing_buffers *p, int wf_size)
{
	int i;
	for(i=0; i<2*p->sample_rate; i++)
		p->waveform[i] = 0;
    
	for(i=0; i < wf_size; i++) {
		float bin[(int)ceil(1 + p->sample_count / wf_size)];
		int j;
		double k = fmod(i+p->phase,wf_size);
		for(j=0;;j++) {
			int n = round(k+j*wf_size);
			if(n >= p->sample_count) break;
			bin[j] = p->samples[n];
		}
		p->waveform[i] = tmean(bin,j);
	}

int step = ceil(wf_size / 100);
	for(i=0; i * step < wf_size; i++)
		p->waveform_sc[i] = p->waveform[i * step];
	quickselect(p->waveform_sc, i, i/2);
	double nl = p->waveform_sc[i/2];
	for(i=0; i<wf_size; i++)
		p->waveform[i] -= nl;

p->waveform_max = vmax(p->waveform, 0, wf_size, &p->waveform_max_i);
}
```

###  Amp 幅值/摆幅计算

```c
/*
平滑窗口
in: 输入数组。
out: 输出数组。
window: 窗口大小。
size: 输入数组的大小。
适合流式数据处理
*/
static void smooth(float *in, float *out, int window, int size)
{
	int i;
	double k = 1 - (1. / window); 	// 衰减因子, 
	double r_av = 0;				// 滑动窗口内元素的累加值。
	double u = 0;					// 当前窗口的最大值估计
	double w = 0;					// 窗口内的最小值估计
    
    // 初始窗口计算, 更新 u 和 r_av
	for(i=0; i < window; i++) {
		u *= k;				// 使用指数衰减来更新窗口内的最大值 u, 减少了对突然变化的响应
		float x = in[i];
		if(x > u) u = x;
		r_av += u;			// 累加器
	}

// 滑动窗口
	for(i=0; i + window < size; i++) {
		out[i] = r_av;
        // 使用指数衰减（通过变量 k）来更新窗口内的最大值 u 和最小值 w，减少重复计算
		u *= k;	  	
		w *= k;
		float x = in[i+window];
		float y = in[i];
		if(x > u) u = x;
		if(y > w) w = y;
		r_av += u - w;
	}
}
```

```c
/*

la  -----  Lift angle 升角：摆轮升角就是摆轮上的圆盘钉与擒纵叉一侧叉口槽壁接触的瞬间，到圆盘钉与擒纵叉另一侧叉口槽壁脱离的瞬间，摆轮所转过的角度。这个角度通常在 48 到 54 度之间。
AMP ----- 摆幅指摆轮的摆动角度，从一个极限位置到另一个极限位置的总角度。用来衡量手表的精度和性能
理想的摆幅通常在 270 到 315 度 之间。这个范围可以确保手表运转正常，具有良好的精度和稳定性。过低或过高的摆幅可能导致误差增大或动力不足。
“基角法” 测量摆幅？？？？
*/
static void compute_amplitude(struct processing_buffers *p, double la)
{
	int i,j,k;
	// 1. 数据初始化处理
	int wf_size = ceil(p->period);
	int window = p->sample_rate / 1000;
	for(i = 0; i < window; i++)
		p->waveform[i + wf_size] = p->waveform[i];
    
    // 2.平滑处理
	float smooth_wf[wf_size];
	smooth(p->waveform, smooth_wf, window, wf_size + window);

// 3.初始化寻找局部最大值，用j记录位置；
	double max = 0;
	for(k = 0; k < 2; k++) {
        // 计算平滑后波形起始位置
		j = floor(fmod((k ? p->tic : p->toc) + p->period/8, p->period));
        
		for(i = 0; i < p->period/8; i++) {
			if(smooth_wf[j] > max) max = smooth_wf[j];
            
			if(++j > p->period) j = 0;
		}
	}
    
    // 4.计算阈值
	double glob_max = vmax(smooth_wf, 0, ceil(p->period), NULL);
	double threshold = fmax(.01 * glob_max, 1.4 * max);  // 阈值
	debug("amp threshold from %s\n", .01 * glob_max > 1.4 * max ? "global maximum" : "noise level");

// 5. 迭代寻找有效振幅，
    // 在阈值小于全局最大值的 20% 时，尝试调整阈值来找到有效脉冲
	p->amp = -1;
	p->tic_pulse = p->toc_pulse = -1;
	while(threshold < .2 * glob_max) {
		debug("amp threshold = %f%% glob max\n", threshold * 100 / glob_max);
		double tic_pulse = -1;
		double toc_pulse = -1;
		for(k = 0; k < 2; k++) {  // 处理两个不同的时间点 tic 和 toc
            // 存储当前段最大值
			double max = 0; 
			// 计算起始位置
			j = floor(fmod((k ? p->tic : p->toc) + 7*p->period/8, p->period)); 
            
            // 查找信号中第一个超过阈值的位置并用j记录位置
			for(i = 0; i < p->period/8; i++) {
				if(smooth_wf[j] > threshold) break;
				if(++j > p->period) j = 0;
			}
            
            // 寻找信号峰值
            // 从超过阈值的点开始，寻找该段的最大值
			for(; i < p->period/8; i++) {
				double x = smooth_wf[j];
				if(x > max) max = x;
				else break;  // 更新最大值，直到信号开始下降
				if(++j > p->period) j = 0;
			}
            
            // 如果找到有效的脉冲段，计算脉冲时间
			if(i < p->period/8) {
                // 计算脉冲时间， period为完整周期的时间 ???
				double pulse = p->period/8 - i - 1;
                
                // 处理两个不同的时间点 tic 和 toc
				if(k) tic_pulse = pulse;
				else toc_pulse = pulse;
				
				debug("amp %s pulse = %f\n", k ? "tic" : "toc", 1000 * pulse / p->sample_rate);
			} else
				goto next_threshold;
		}
        
        // 6. 计算振幅和脉冲
		double tic_amp_abs = .5 / sin(M_PI * tic_pulse / p->period);
		double toc_amp_abs = .5 / sin(M_PI * toc_pulse / p->period);
        // 6.1 用升角 la 修正计算结果
		double tic_amp = la * tic_amp_abs;  
		double toc_amp = la * toc_amp_abs;
        
        // 7.判断振幅是否在合理范围内，并且两个振幅的差异要求小于60
		if(135 < tic_amp && tic_amp < 360 && 135 < toc_amp && toc_amp < 360 && fabs(tic_amp - toc_amp) < 60) 		{
			p->amp = (tic_amp_abs + toc_amp_abs) / 2;
			p->tic_pulse = tic_pulse;
			p->toc_pulse = toc_pulse;
            // be(偏振):计算摆轮在周期中相对于中心点的偏差，用于衡量摆轮的偏移程度
			p->be = p->period/2 - fabs(p->toc - p->tic + p->tic_pulse - p->toc_pulse);
			debug("amp: be = %.1f\n",fabs(p->be)*1000/p->sample_rate);
			debug("amp = %f\n", la * p->amp);
			break;
		} else
			debug("amp rejected\n");

// 8. 重新调整阈值
            next_threshold:	
                threshold *= 1.4;
		}
	if(p->amp < 0) debug("amp failed, amp threshold = %f, .2 * glob_max = %f\n", threshold, .2 * glob_max);
}
```

# 参考文档

[1. 节拍误差的定义](https://www.watchuseek.com/threads/definition-of-beat-error.2394130/)

[2. A program for timing mechanical watches doc](https://www.watchuseek.com/threads/open-source-timing-software.2542874/)

[3. 机械表的工作原理](https://www.watchuseek.com/threads/how-a-mechanical-watch-works.772681/)

[4. 擒纵机构工作原理解析](https://www.bilibili.com/video/BV1z44y1o734/?spm_id_from=333.337.search-card.all.click&vd_source=64bcdab1e8263295a838d5e4a0fab29c)

[5.摆轮升角](https://baike.xbiao.com/291.html)

[6.机械手表校表仪的打点线有什么含义](https://www.bilibili.com/video/BV1EK411C7i6/?spm_id_from=333.337.search-card.all.click&vd_source=64bcdab1e8263295a838d5e4a0fab29c)

[7. 机械表滴答滴答声音的含义](https://www.bilibili.com/video/BV16d4y157xQ/?spm_id_from=333.788.recommend_more_video.-1&vd_source=64bcdab1e8263295a838d5e4a0fab29c)

[擒纵机构：机械之美](https://www.sic.cas.cn/zt/kpwsx/kpwz/201509/t20150907_4421143.html)

[**适用于 Windows 的 7 款最佳免费音频频谱分析仪软件**](https://www.audioapp.cn/bbs/thread-219136-1-1.html)