使用格兹尔算法的循环分析(基础篇)
◍ 用格兹尔算法拆价格周期
格兹尔(Goertzel)算法原本用于数字信号处理里提取特定频率分量,放到 MT5 里可以拿来做循环分析,把看似随机的汇率波动拆出隐藏的周期结构。这套思路在 2024 年 2 月的一篇指标实现里被搬进 MQL5,发布后约一个月获得 1085 次查看、3 条反馈,说明手动交易者对这个轻量频域工具确有需求。 和傅里叶变换一次性算全频谱不同,格兹尔只盯你指定的几个频率,计算量小、适合实时跑在图表上。对外汇和贵金属这类高杠杆品种,周期识别只提供概率倾向,实际进出场仍要叠加价格行为确认,杠杆风险不可轻视。 想验证的话,直接在 MT5 新建指标工程,把下方核心循环贴进去,改一下目标频率参数就能看到不同周期分量的强度变化。
class="type">class="kw">double Goertzel(class="type">class="kw">double &data[], class="type">int N, class="type">class="kw">double freq) { class="type">class="kw">double omega = class="num">2.0 * M_PI * freq / N; class="type">class="kw">double coeff = class="num">2.0 * cos(omega); class="type">class="kw">double s0 = class="num">0.0, s1 = class="num">0.0, s2 = class="num">0.0; for(class="type">int i = class="num">0; i < N; i++) { s0 = data[i] + coeff * s1 - s2; s2 = s1; s1 = s0; } class="kw">return s1 * s1 + s2 * s2 - coeff * s1 * s2; }
格兹尔算法能挖出价格里的主周期
格兹尔算法(Goertzel algorithm)本质是一种数字信号处理手段,专长于在信号里抠出某一个特定频率分量,比起做全套 FFT 要省算力和延迟。它精度在线、能实时跑,这两点正好对上金融时间序列的脾气——K线说到底也是带噪声的离散序列。 这篇文章里我们要做的,是用它从报价流里反推出当前行情的主导周期。主导周期不是玄学,是价格波动里能量最集中的那个频段,抓到它,后续做均值回归或周期跟随策略才有锚点。 后面会直接给 MQL5 的实现代码,并演示怎么调用去识别价格报价中的周期成分。你打开 MT5 自建一个指标把代码贴进去,就能看到不同窗口长度下算出的主频偏移,外汇和贵金属杠杆高、跳空多,周期容易被打断,验证时建议先用 XAUUSD 的 M5 跑一遍看稳定性。
「用格兹尔算法抠出价格里的特定频率」
格兹尔算法由 Gerald Goertzel 在 1958 年提出,本质是用递推方式算离散傅立叶变换(DFT)的单项。和 FFT 比,它只在盯少数几个频率分量时更省算力——这对行情数据很实用,因为我们往往只关心某条周期有没有能量。 公式层面,它在频率仓 k 上的累积幅度靠两个历史状态往前滚:X[k]=X[k-1]+x[n]-coeff*X[k-2],其中 coeff=2*cos(2πk/N),N 是样本数也就是看的柱数。迭代完所有样本后,最后三个状态值直接给出该频率 bin 的实部与虚部,不需要全套 DFT。 DFT 能分辨的频率下限是 1/N、上限是 (N/2)/N,步长卡死在 1/N。也就是说看 200 根柱,最小可探周期就是 200 根,中间频带探不到——这正是 Ehlers 搞 MESA 的初衷。 但 D.Meyers 的论文指出,噪声大的金融序列里格兹尔反而可能比 MESA 谱分辨更好。外汇和贵金属报价天然高噪,拿这套做周期检测倾向更稳,但杠杆品种高风险,参数错了就是瞎猜。
◍ 用 Goertzel 只抓你关心的周期带
Goertzel 算法不像全段 FFT 那样把整个频谱都算出来,它一次只对有限个频率做采样,所以吃 CPU 少、出数快。做价格周期分析时,我们往往只盯 8~50 根 K 线的波段,没必要算全频。 这个类靠最小/最大周期来框定频带:周期越小频率越高,周期越大频率越低。你可以用 SetMinMaxPeriodLength() 提前设好带宽,也可以在调 Dft() 时顺手传 min_wavePeriod 和 max_wavePeriod 覆盖。 Dft 有四个重载:两种输出复数阵列(complex[]),两种把实部、虚部分开吐到 double[]。实虚部分离的版本在写自定义振幅/相位指标时更好接,不用再拆结构体。 注意 SetMinMaxPeriodLength 里有硬校验:min_wavePeriod 不能小于 2,也不能大于等于 max_wavePeriod,否则直接 Print 报错并返回 false。外汇和贵金属波动里周期带设错,出来的分量就是噪声,杠杆品种高风险,参数先拿历史数据验证再上实盘。 下面这段是类声明和设周期函数的骨架,逐行拆一下关键行: // class CGoertzel:封装 Goertzel 算法,private 里 m_minPeriod/m_maxPeriod 存频带边界 // 构造函数 CGoertzel(void):把上下周期都置 0,等后面显式设 // GetMinPeriodLength / GetMaxPerodLength:读当前频带(注意原拼写 MaxPerod 少了 i) // SetMinMaxPeriodLength:校验 min<2 或 min>=max 就打印错误并返回 false,否则写成员返回 true // Dft 四个重载:前两个用类内已设频带,后两个调用时传频带;输出可选 complex 或 实/虚双数组
class=class="str">"cmt">//+------------------------------------------------------------------+ class=class="str">"cmt">//| CGoertze class implementing Goertzel algorithm | class=class="str">"cmt">//+------------------------------------------------------------------+ class CGoertzel { class="kw">private: class="type">uint m_minPeriod; class="type">uint m_maxPeriod; class="kw">public: class=class="str">"cmt">//constructor CGoertzel(class="type">void) { m_minPeriod=m_maxPeriod=class="num">0; } class=class="str">"cmt">//destructor ~CGoertzel(class="type">void) { } class=class="str">"cmt">//get methods class="type">uint GetMinPeriodLength(class="type">void) { class="kw">return m_minPeriod;} class="type">uint GetMaxPerodLength(class="type">void) { class="kw">return m_maxPeriod;} class=class="str">"cmt">//set methods class="type">bool SetMinMaxPeriodLength(const class="type">uint min_wavePeriod, const class="type">uint max_wavePeriod); class=class="str">"cmt">// goertzel dft transform methods class="type">bool Dft(const class="type">class="kw">double &in_data[],complex &out[]); class="type">bool Dft(const class="type">class="kw">double &in_data[], class="type">class="kw">double &out_real[], class="type">class="kw">double &out_imaginary[]); class="type">bool Dft(const class="type">uint min_wavePeriod,const class="type">uint max_wavePeriod,const class="type">class="kw">double &in_data[],complex &out[]); class="type">bool Dft(const class="type">uint min_wavePeriod,const class="type">uint max_wavePeriod,const class="type">class="kw">double &in_data[], class="type">class="kw">double &out_real[], class="type">class="kw">double &out_imaginary[]); }; class=class="str">"cmt">//+------------------------------------------------------------------+ class=class="str">"cmt">//| Set the Minimum and maximum periods of selected frequency band | class=class="str">"cmt">//+------------------------------------------------------------------+ class="type">bool CGoertzel::SetMinMaxPeriodLength(const class="type">uint min_wavePeriod,const class="type">uint max_wavePeriod) { if(min_wavePeriod<class="num">2 || min_wavePeriod>=max_wavePeriod) { Print("Critical error min_wavePeriod cannot be less than max_wavePeriod or less than class="num">2"); class="kw">return false; } m_minPeriod=min_wavePeriod; m_maxPeriod=max_wavePeriod; class="kw">return true; } class=class="str">"cmt">//+-----------------------------------------------------------------------------------+ class=class="str">"cmt">//| Calculate the frequency domain representation of input array as complex array | class=class="str">"cmt">//+-----------------------------------------------------------------------------------+
Goertzel 类的 DFT 重载与频域计算内核
CGoertzel 类提供了三组 Dft 重载入口,前两个只是带周期范围参数的包装函数,内部都先调用 SetMinMaxPeriodLength 做边界设定,再转交核心实现。真正干活的第三个重载把输入序列转成复数数组,输出每个周期分量的实部与虚部。 核心函数一上来就卡样本量:要求输入数组长度至少为最大周期值的 3 倍(minsize = 3 * m_maxPeriod),否则直接 Print 报错并返回 false。同时强制 m_minPeriod 必须 ≥2 且小于 m_maxPeriod,这两个约束不合规也会中断,避免后续除零或无效频点。 逐周期扫描时,索引 i 小于 m_minPeriod 的分量直接清零跳过;其余分量用 Goertzel 递推式 v0 = coeff*v1 - v2 + in_data[k] 从后往前滚一遍,coeff 由 2*cos(2π/i) 算出。最终实部取 v1 - v2*0.5*coeff,虚部取 v2*sin(2π/i),写进 out[i]。 在 MT5 里把这段直接塞进 EA 或脚本,接一段 EURUSD 的收盘价数组跑一下,就能看到不同周期(比如 10~50 根 K 线)对应的频域幅值。外汇与贵金属波动受杠杆与事件驱动影响,频域峰值只代表历史周期倾向,实盘使用前务必用小资金验证高风险敞口。
class="type">bool CGoertzel::Dft(const class="type">uint min_wavePeriod,const class="type">uint max_wavePeriod,const class="type">class="kw">double &in_data[],complex &out[]) { if(!SetMinMaxPeriodLength(min_wavePeriod,max_wavePeriod)) class="kw">return(false); class="kw">return Dft(in_data,out); } class=class="str">"cmt">//+------------------------------------------------------------------------------------+ class=class="str">"cmt">//| Calculate the frequency domain representation of input array as two separate arrays| class=class="str">"cmt">//+------------------------------------------------------------------------------------+ class="type">bool CGoertzel::Dft(const class="type">uint min_wavePeriod,const class="type">uint max_wavePeriod,const class="type">class="kw">double &in_data[], class="type">class="kw">double &out_real[], class="type">class="kw">double &out_imaginary[]) { if(!SetMinMaxPeriodLength(min_wavePeriod,max_wavePeriod)) class="kw">return(false); class="kw">return Dft(in_data,out_real,out_imaginary); } class=class="str">"cmt">//+-----------------------------------------------------------------------------------+ class=class="str">"cmt">//| Calculate the frequency domain representation of input array as complex array | class=class="str">"cmt">//+-----------------------------------------------------------------------------------+ class="type">bool CGoertzel::Dft(const class="type">class="kw">double &in_data[],complex &out[]) { class="type">uint minsize=(class="num">3*m_maxPeriod); class="type">uint fullsize=in_data.Size(); if(fullsize<minsize) { Print("Sample size too small in relation to the largest period cycle parameter"); class="kw">return false; } if(m_minPeriod>=m_maxPeriod || m_minPeriod<class="num">2) { Print("Critical error: Invalid input parameters :- max_period should be larger than min_period and min_period cannot be less than class="num">2"); class="kw">return false; } if(out.Size()!=m_maxPeriod) ArrayResize(out,m_maxPeriod); class="type">class="kw">double v0,v1,v2,freq,coeff,real,imag; for(class="type">uint i=class="num">0; i<m_maxPeriod; i++) { if(i<m_minPeriod) { out[i].imag=out[i].real=class="num">0.0; class="kw">continue; } v0=v1=v2=class="num">0.0; freq=MathPow(i,-class="num">1); coeff=class="num">2.0*MathCos(class="num">2.0*M_PI*freq); for(class="type">uint k=minsize-class="num">1; k>class="num">0; k--) { v0=coeff*v1-v2+in_data[k]; v2=v1; v1=v0; } real=v1-v2*class="num">0.5*coeff; imag=v2*MathSin(class="num">2*M_PI*freq); out[i].real=real; out[i].imag=imag; } class="kw">return true; } class=class="str">"cmt">//+------------------------------------------------------------------------------------+ class=class="str">"cmt">//| Calculate the frequency domain representation of input array as two separate arrays| class=class="str">"cmt">//+------------------------------------------------------------------------------------+
「Goertzel 离散傅里叶变换的边界与递推实现」
在 MT5 里跑周期检测,样本长度先卡死一条线:最小样本量必须是最大周期参数的 3 倍(minsize = 3 * m_maxPeriod)。若传入序列长度不足,函数直接 Print 报错并返回 false,不会吐出任何频谱——这是很多自定义指标在周线级品种上静默失效的根因。 min_period 与 max_period 的合法性也绕不开:max 必须大于 min,且最小周期不能低于 2。否则同样走 false 分支。这两道校验放在循环前,等于把脏参数挡在算力消耗之外。 核心递推在双循环里。外层按 i 从 0 扫到 m_maxPeriod,小于 min_period 的频点直接置零跳过;大于等于时才进入 Goertzel 单频点算法。内层从 minsize-1 倒序累加到 1,用 v0 = coeff*v1 - v2 + in_data[k] 做二阶差分递推,coeff 由 2*cos(2π/i) 给出。 实部和虚部的提取只用末尾两个状态量:real = v1 - v2*0.5*coeff,imag = v2*sin(2π/i)。比起完整 FFT,这套写法在只盯某几个候选周期时更省,但外汇与贵金属波动有跳空噪声,周期识别结果仅作概率参考,实盘须自担高风险。
class="type">bool CGoertzel::Dft(const class="type">class="kw">double &in_data[],class="type">class="kw">double &out_real[],class="type">class="kw">double &out_imaginary[]) { class="type">uint minsize=(class="num">3*m_maxPeriod); class="type">uint fullsize=in_data.Size(); if(fullsize<minsize) { Print("Sample size too small in relation to the largest period cycle parameter"); class="kw">return false; } if(m_minPeriod>=m_maxPeriod || m_minPeriod<class="num">2) { Print("Critical error: Invalid input parameters :- max_period should be larger than min_period and min_period cannot be less than class="num">2"); class="kw">return false; } if(out_real.Size()!=m_maxPeriod) ArrayResize(out_real,m_maxPeriod); if(out_imaginary.Size()!=m_maxPeriod) ArrayResize(out_imaginary,m_maxPeriod); class="type">class="kw">double v0,v1,v2,freq,coeff,real,imag; for(class="type">uint i=class="num">0; i<m_maxPeriod; i++) { if(i<m_minPeriod) { out_real[i]=out_imaginary[i]=class="num">0.0; class="kw">continue; } v0=v1=v2=class="num">0.0; freq=MathPow(i,-class="num">1); coeff=class="num">2.0*MathCos(class="num">2.0*M_PI*freq); for(class="type">uint k=minsize-class="num">1; k>class="num">0; k--) { v0=coeff*v1-v2+in_data[k]; v2=v1; v1=v0; } real=v1-v2*class="num">0.5*coeff; imag=v2*MathSin(class="num">2*M_PI*freq); out_real[i]=real; out_imaginary[i]=imag; } class="kw">return true; }