|
|
发表于 2025-6-16 18:33:42
|
显示全部楼层
本帖最后由 snakeemail 于 2025-6-16 18:36 编辑
我把相位改到90°，计算的相位phase: 359.686 degree
#include <math.h>#include <complex.h>#include <stdio.h>
void optimized_goertzel(double *x, int N, double k, double *magnitude, double *phase) { double w = 2.0 * M_PI * k / N; double cosine = cos(w); double sine = sin(w); double coeff = 2.0 * cosine;
double q0, q1 = 0.0, q2 = 0.0;
for (int n = 0; n < N; n++) { q0 = coeff * q1 - q2 + x[n]; q2 = q1; q1 = q0; }
// 直接计算实部和虚部 double real = q1 - q2 * cosine; double imag = q2 * sine;
*magnitude = sqrt(real*real + imag*imag); *phase = atan2f(imag, real); if (imag < 0 ) { *phase = 360 + *phase; }}
void test_goertzel() {
// 示例信号参数 const int N = 1000;//200; // 采样点数 const double fs = 1000; // 采样率(Hz) const double target_freq = 50; // 目标频率(Hz)
// 生成测试信号: 50Hz正弦波 + 少量噪声 double signal[N]; for (int i = 0; i < N; i++) { double t = i / fs; signal[i] = sin(2 * M_PI * 50 * t + (90 * M_PI / 180)); // + 0.1 * sin(2 * M_PI * 120 * t); // }
// 计算对应的DFT bin index double k = (target_freq * N / fs); //50 * 1000 / 1000
// 使用Goertzel算法计算 //double complex result = goertzel(signal, N, k); double mag_peak = 0; // = compute_magnitude(result); double phase = 0; // = compute_phase(result);
optimized_goertzel(signal, N, k, &mag_peak, &phase); // 输出结果 printf("目标频率: %.1f Hz\n", target_freq); printf("对应DFT bin: %.2f\n", k); //printf("复数结果: %.3f + %.3fi\n", creal(result), cimag(result)); printf("amp: %.3f\n", mag_peak / (N/2)); printf("phase: %.3f degree\n", phase);}
int main(void) //int argc, char *argv[]){ test_goertzel();
return 0;} |
|
楼主
。 
