hobbylab.jp Software algorithm Goertzel

ゲーツェルアルゴリズムは、一言で言うと「特定の周波数だけに激しく共振する『デジタル上の音叉（おんさ）』を作る技術」です。
FFT（高速フーリエ変換）が全体の周波数をまとめて計算するのに対し、ゲーツェルは「700Hzだけ」を狙い撃ちするため、マイコンにとても優しい（計算が軽い）のが特徴です。
無線通信のモールス信号のディコーダを作成するために使用したアルゴリズムです。
したがってここで説明するものはモールス信号のディコーダを作成するための説明となります。

(1)モールス信号とは

・基本の2種類: 短い「トン（・）」と、その3倍の長さの長い「ツー（－）」を使います。
・間隔のルール: 符号の間は短点1つ分、文字の間は短点3つ分、単語の間は短点7つ分あけます。
・有名な例: 遭難信号の「SOS」は「・資料・－・－・（トトト・ツーツー・トトト）」という分かりやすい形です。
波形で示すと

となりますが、この信号を音で聞き取るため、黄色の部分に600〜800(hz)の音を付けています。
以下からはこの周波数を700(Hz)と仮定して説明していきます。
信号 S の部分を下に示します。(紫色の交流波形が例えば700(Hz)のつもりです!)

(2)ゲーツェルアルゴリズムはなにをするか

ゲーツェルアルゴリズムは(1)で説明した700(Hz)の信号を(音叉のように)限定して信号の外枠を取り出すために使用します。

(3) 概略装置

どんな装置で構成しているか説明します。

①ラジオー> ②アンプ・フィルター> ③ MCU ー> ④ LED (オシロで確認)
※将来は
④フィルター> ⑤アンプー> ⑥ミニブザーでモニターできる回路にしたい
及び
④表示器

①ラジオは Si4735 を使用します。
②アンプ・フィルターはオペアンプ TVL9062 を使用します。
③MCU は CH32V203C8T6 を使用します。
④今回はオシロスコープで整形した波形を確認します。
将来
④フィルターは R (抵抗) C (コンでサー)で構成したフィルターです。
⑤アンプはLTK5128を使用します。
⑥ミニブザーはパッシブ型を使用します。
および
④表示器はLCDでモールス符号を表示したいと思います。

(4)ゲーツェルアルゴリズムのプログラム

今回説明するのは ③ の部分のプログラムになります。
プログラムは大きく分けて3つのステップでモールス信号の音を限定しています。

【入力信号】(PB0からADC)
　　│
　　▼
［ステップ1: 直流カット］ ─── 処理場所: x = adc_buffer[i] - CentralSignal;
　　│ (1.533Vの直流成分を引き算して、純粋な交流の波にする)
　　▼
［ステップ2: 音叉の共振］ ─── 処理場所: s = x + coeff * s_prev - s_prev2;
　　│ (TargetFreqで設定した周波数が来ると、s の値がどんどん巨大化する)
　　▼
［ステップ3: パワー測定］ ─── 処理場所: power = (s_prev^2) + (s_prev2^2) - ...
　　│ (溜まった振動のエネルギーを計算し、ExclusionSignalより大きければPB8をHIGHに)
　　▼
【結果出力】(PB8ピン)

1.2 動作状況

モールス信号は ③MCU の PB0 端子に接続されます。
また並列にオシロスコープに接続します。(下図の水色線)
ゲーツェルアルゴリズムを使用したプログラムによりモールス信号の外枠を取り出します。
その信号が ③MCU の PB8 から出力され ④オシロスコープに接続されています。(下図の黄色線)

2.プログラム

2.1 内容

周波数は前段のフィルター特性から690(Hz)とした。
サンプリングは約10倍の7000(Hz)とした。

#include "debug.h"
#include <math.h>

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

// 定数定義
#define SAMPLING_RATE 7000.0f  // サンプリング周波数 (6~8kHzの中央)
#define BLOCK_SIZE    140      // ブロックサイズ (140サンプル = 20msフレーム)

// ユーザー指定の制御変数
volatile float TargetFreq = 690.0f;     // 690Hzに固定（実行中に600Hz〜800Hzで調整も可能）
volatile int16_t CentralSignal = 2026;  // 12bit ADCの無信号時中心値 (交流0点=1.533V.)
volatile float ExclusionSignal = 110000.0f; //

// ゲーツェルアルゴリズム計算用ワーク変数
volatile int sample_count = 0;
float adc_buffer[BLOCK_SIZE];
volatile uint8_t processing_flag = 0;
volatile int32_t power_count = 0;   // カウンタ
volatile int64_t apower_sum = 0;    // 200回分の合計を蓄える巨大な変数

// 周辺機能の初期化
void GPIO_Config(void) {
    GPIO_InitTypeDef GPIO_InitStructure = {0};

// クロック有効化 (GPIOB)
    RCC_APB2PeriphClockCmd(RCC_APB2Periph_GPIOB, ENABLE);

// PB0: アナログ入力 (ADC1_IN8)
    GPIO_InitStructure.GPIO_Pin = GPIO_Pin_0;
    GPIO_InitStructure.GPIO_Mode = GPIO_Mode_AIN;
    GPIO_Init(GPIOB, &GPIO_InitStructure);

// PB8: デジタル出力 (検出フラグ用)
    GPIO_InitStructure.GPIO_Pin = GPIO_Pin_8;
    GPIO_InitStructure.GPIO_Mode = GPIO_Mode_Out_PP; // プッシュプル出力
    GPIO_InitStructure.GPIO_Speed = GPIO_Speed_50MHz;
    GPIO_Init(GPIOB, &GPIO_InitStructure);

GPIO_WriteBit(GPIOB, GPIO_Pin_8, Bit_RESET);
}

void ADC_Config(void) {
    ADC_InitTypeDef ADC_InitStructure = {0};

// クロック有効化 (ADC1)
    RCC_APB2PeriphClockCmd(RCC_APB2Periph_ADC1, ENABLE);
    RCC_ADCCLKConfig(RCC_PCLK2_Div6); // ADCクロック設定

ADC_DeInit(ADC1);
    ADC_InitStructure.ADC_Mode = ADC_Mode_Independent;
    ADC_InitStructure.ADC_ScanConvMode = DISABLE;
    ADC_InitStructure.ADC_ContinuousConvMode = DISABLE; // タイマーでトリガーするため単発モード
    ADC_InitStructure.ADC_ExternalTrigConv = ADC_ExternalTrigConv_None;
    ADC_InitStructure.ADC_DataAlign = ADC_DataAlign_Right;
    ADC_InitStructure.ADC_NbrOfChannel = 1;
    ADC_Init(ADC1, &ADC_InitStructure);

ADC_Cmd(ADC1, ENABLE);

// 校正
    ADC_ResetCalibration(ADC1);
    while(ADC_GetResetCalibrationStatus(ADC1));
    ADC_StartCalibration(ADC1);
    while(ADC_GetCalibrationStatus(ADC1));
}

// PB0(Channel 8)からA/D変換値を取得する関数
uint16_t Get_ADC_Val(void) {
    ADC_RegularChannelConfig(ADC1, ADC_Channel_8, 1, ADC_SampleTime_239Cycles5);
    ADC_SoftwareStartConvCmd(ADC1, ENABLE);
    while(!ADC_GetFlagStatus(ADC1, ADC_FLAG_EOC));
    return ADC_GetConversionValue(ADC1);
}

// 8kHzサンプリング用タイマー設定 (TIM2)
void TIM2_Config(void) {
    TIM_TimeBaseInitTypeDef TIM_TimeBaseStructure = {0};
    NVIC_InitTypeDef NVIC_InitStructure = {0};

RCC_APB1PeriphClockCmd(RCC_APB1Periph_TIM2, ENABLE);

// 144MHz動作を想定時の7kHz設定 (144,000,000 / 144 / 144 = 6944Hz)
    TIM_TimeBaseStructure.TIM_Period = 144 - 1;
    TIM_TimeBaseStructure.TIM_Prescaler = 144 - 1;

TIM_TimeBaseStructure.TIM_ClockDivision = TIM_CKD_DIV1;
    TIM_TimeBaseStructure.TIM_CounterMode = TIM_CounterMode_Up;
    TIM_TimeBaseInit(TIM2, &TIM_TimeBaseStructure);

TIM_ITConfig(TIM2, TIM_IT_Update, ENABLE);

NVIC_InitStructure.NVIC_IRQChannel = TIM2_IRQn;
    NVIC_InitStructure.NVIC_IRQChannelPreemptionPriority = 1;
    NVIC_InitStructure.NVIC_IRQChannelSubPriority = 1;
    NVIC_InitStructure.NVIC_IRQChannelCmd = ENABLE;
    NVIC_Init(&NVIC_InitStructure);

TIM_Cmd(TIM2, ENABLE);
}

// Goertzelフィルター本体の計算関数
uint8_t Process_Goertzel(void) {
    float s_prev = 0.0f;
    float s_prev2 = 0.0f;
    float s;

// 現在の TargetFreq から係数を計算
    float omega = (2.0f * M_PI * TargetFreq) / SAMPLING_RATE;
    float coeff = 2.0f * cosf(omega);

// IIR フィルタの処理
    for (int i = 0; i < BLOCK_SIZE; i++) {
        // 入力からCentralSignal(直流オフセット)を引いて純粋な交流成分にする
        float x = adc_buffer[i] - (float)CentralSignal;

s = x + coeff * s_prev - s_prev2;
        s_prev2 = s_prev;
        s_prev = s;
    }

// パワー（振幅の2乗）の計算
    float power = (s_prev * s_prev) + (s_prev2 * s_prev2) - (coeff * s_prev * s_prev2);

// 計算の正規化（必要に応じて調整）
    power = power / (BLOCK_SIZE * BLOCK_SIZE);

// ─── 【修正版】power モニター ───
    power_count++;
    apower_sum += (int32_t)power;

if (power_count >= 200){
        // 合計値を回数で割り算して平均を出す（%ld で安全に出力）
        printf("After power (Avg) = %ld \r\n", (long)(apower_sum / power_count));
        power_count = 0;
        apower_sum = 0;
    }
    // ─── power モニター END ───

// 雑音閾値(ExclusionSignal)を超えているか判定
    if (power > ExclusionSignal) {
        return 1; // 検出成功
    } else {
        return 0; // 不検出 (雑音以下)
    }
}

/*********************************************************************
 * @fn      main
 *
 * @brief   Main program.
 *
 * @return  none
 */
int main(void)
{
    NVIC_PriorityGroupConfig(NVIC_PriorityGroup_1);
    SystemCoreClockUpdate();
    Delay_Init();
    USART_Printf_Init(115200);
    printf("SystemClk:%d\r\n", SystemCoreClock);
    printf( "ChipID:%08x\r\n", DBGMCU_GetCHIPID() );
    printf("This is printf example\r\n");

GPIO_Config();
    ADC_Config();
    TIM2_Config();

while(1) {
        // バッファが一杯になったらメインループでゲーツェル演算を行う
        if (processing_flag) {
            uint8_t detected = Process_Goertzel();

if (detected) {
                // 信号を検出した場合はPB8をHIGH
                GPIO_WriteBit(GPIOB, GPIO_Pin_8, Bit_SET);
            } else {
                // 検出されなかった場合はPB8をLOW
                GPIO_WriteBit(GPIOB, GPIO_Pin_8, Bit_RESET);
            }

// フラグを戻してサンプリングを再開
            processing_flag = 0;
        }
    }
}

// タイマー割り込みハンドラ (8kHz周期で呼ばれる)
void TIM2_IRQHandler(void) __attribute__((interrupt("WCH-Interrupt-fast")));
void TIM2_IRQHandler(void) {
    if (TIM_GetITStatus(TIM2, TIM_IT_Update) != RESET) {
        TIM_ClearITPendingBit(TIM2, TIM_IT_Update);

// バッファ処理中でなければサンプリング
        if (!processing_flag) {
            adc_buffer[sample_count] = (float)Get_ADC_Val();
            sample_count++;

if (sample_count >= BLOCK_SIZE) {
                sample_count = 0;
                processing_flag = 1; // メイン側に計算要求を出してサンプリングを一時待機
            }
        }
    }
}

2.2 説明

プログラムが中で何をしているのかを「ブランコと音叉（おんさ）」のイメージで、コードの1行1行と照らし合わせながら直感的に説明します。

ゲーツェル関数の全体イメージ：デジタルの中の「ブランコ」
この関数は、入ってきた波（音声信号）のタイミングに合わせて、デジタル空間にある「ブランコ」をタイミングよく押し続ける処理をしています。
●もし入力された音が「狙った周波数（700Hz）」なら、ブランコを押すタイミングがピッタリ合うので、ブランコはどんどん高く揺れます（数値が巨大化する）。
●もし違う周波数の音やノイズなら、押すタイミングがバラバラで打ち消し合うため、ブランコはほとんど揺れません（数値が大きくならない）。

このイメージを持って、実際のコード「Process_Goertzel( ) 」を見ていきましょう。

2.2.1 変数１

●何をしているか: ブランコの「過去の揺れ具合」を記録するメモ帳を準備しています。
●仕組み: ゲーツェルでは、「1回前の揺れ（s_prev）」と「2回前の揺れ（s_prev2）」の2つだけを記憶しながら計算を進めます。
最初は動いていないので 0 を入れます。

2.2.2 変数２

●何をしているか: 狙う周波数（700Hz）にピッタリ合う「ブランコの長さ（係数 coeff）」を計算しています。
●仕組み: ここは数学の数式ですが、マイコンに対して「今から700Hzのリズムを基準にするよ」と教えている定数（決まった数字）を作っているだけ、と捉えてください。

2.2.3 直流分除外

●何をしているか: 160個溜まったA/D変換のデータ（波の形）を、1個ずつ順番に取り出して「直流（下駄）を脱がせる」処理をしています。
●仕組み: A/D変換の値は常にプラスのボルト（1.533Vなど）を中心に上下しています。
そのまま足し算すると合計値がどんどんプラスに偏ってしまうため、中心値（CentralSignal）を引き算して、純粋な「＋とーの綺麗な交流の波（x）」に変換しています。

2.2.4 ゲーツェルアルゴリズム

●何をしているか: ここがゲーツェルアルゴリズムの心臓部（ブランコを押す瞬間）です！
●仕組み:

・現在の波の高さ（x）に、「1回前の揺れ（s_prev）にリズム係数を掛けたもの」を足し算し、さらに「2回前の揺れ（s_prev2）」を引き算します。
・この「過去のデータを絶妙なタイミングで足し引きする」という算数マジックにより、690Hzの波が来たときだけ、計算結果の s が雪だるま式にどんどん巨大化していきます。
・次の計算のために、メモ帳（s_prev と s_prev2）の内容を最新版に更新して、160回これを繰り返します。

2.2.5 外枠のHi判定

●何をしているか: 160回押し終わった後、最終的に「ブランコがどれだけ激しく揺れているか（エネルギーの強さ）」を計算しています。
●仕組み: 最後の最後で、残ったメモ帳のデータ（1回前と2回前の値）を掛け合わせることで、「700Hzの成分がどれだけ強く含まれていたか」という1つの正の数字（power）を弾き出します。

2.2.6 スケール調整

何をしているか: 数字が大きくなりすぎているので、サンプリング数で割り算して「見慣れたスケール（振幅の2乗）」にサイズを縮小しています

2.2.7 まとめ

ゲーツェルアルゴリズムの中身は、「過去の計算結果をちょっと加工して、次の計算に使い回すループ処理」です。
難しい掛け算や引き算に見えますが、やっていることは「690Hzのリズムに合わせて数字を雪だるま式に大きくする仕掛け」と、最後に「その雪だるまの大きさをチェックして合否を決める仕掛け」の2つだけです。
だからFFTのような複雑な処理をしなくても、マイコンで超高速に動かすことができます。

2.3 補足説明

2.3.1 入力電圧の設定

　(1)センター電圧　CentralSignal = 2026;

A/D変換が12bitで電源電圧（基準電圧）が3.1Vのとき、1bit（最小分解能）は約0.000757V（0.757mV）になります。
1bitあたりの電圧＝ 3.1V / 4096 =0.0007598
入力電圧　中心電圧1.533V時のbit　＝　1.533 × 4096 / 3.1 = 2026

　(2)ノイス除外電圧　ExclusionSignal = 110000.0f;

ノイズ(A/D値) = 0.3 × 4096 / 3.1 = 296.4
しきい値(パワー) = 396.4×396.4 / 2 ＝ 78,400 ≒ 110000

しきい値は「最大ノイズの1.5倍〜2倍」に設定する
安全マージン（余裕）を持たせるため、見つかった最大ノイズパワーの 1.5倍〜2倍を最初のしきい値（ExclusionSignal）に設定します。
(1)上記の例（最大ノイズが 82000）であれば：82000 × 1.5 ＝ 123000 あたりを狙い目にします。
(2)これにより、ノイズが少し突発的に大きくなっても、PB8が誤ってHIGHになるのを防げます。

2.3.2 プログラムの具体的役割

(1)狙う周波数の決定 (coeff = 2.0f * cosf(omega);)

役割: 音叉の「長さ」や「硬さ」を決める工程です。
TargetFreq（700Hz）に合わせて、どの周波数に共振するかをここで数式化しています。

(2)共振ループ (for (int i = 0; i < BLOCK_SIZE; i++))

役割: 入ってくる波（x）に対して、過去のデータ（s_prev, s_prev2）を絶妙なタイミングで足し引きします。
もし入力信号が700Hzなら、ブランコをタイミングよく押すように内部の数値（s）がどんどん増幅していきます。
違う周波数なら打ち消し合って大きくなりません。

(3)判定と出力 (if (power > ExclusionSignal))

役割: ループ終了後、どれだけ激しく共振したかを計算します。
無信号時のノイズによる共振（約78,400）を超えて、本物の信号が来たと確信できた（110,000を超えた）場合のみ、PB8ピンをHIGHにします。

2.3.3 検出時間とブロックサイズ

a.検出を短くする方法

現在の設定（サンプリングレート 7000Hz、BLOCK_SIZE = 140）では、160個のデータが溜まるのを待つため、計算を開始するまでにどうしても 140 ÷ 7000 = 0.02秒 (20ms) の待ち時間が発生します。
例えば、BLOCK_SIZE を 70 に半減またはサンプリングを14000(Hz)に増加させれば、データ収集待ちは 10ms に短縮されます。

b.ブロックサイズは概略説明のどこの部分でどんな処理か。

概略説明の中の「ステップ2: 音叉の共振（forループ）」の回数そのものです。
プログラムは、指定された BLOCK_SIZE 回だけ波を観察し、音叉を繰り返し叩いて「どれだけ大きく揺れたか」をテストしています。
●小さくするメリット: 反応速度が速くなります（遅延が減る）。
●小さくするデメリット: 音叉を叩く回数が減るため、十分に共振せず、周波数の見分け性能（キレ）が悪くなります。また、ノイズの影響を受けやすくなります。

2.3.4 共振周波数の変更

コントロールするのがまさに「ブロックサイズ（BLOCK_SIZE）」です。

ゲーツェルアルゴリズムには「帯域幅（Q値）を直接可変するパラメータ」はありません。
その代わり、ブロックサイズを大きくするほど周波数の選択性がシビア（ピンポイント）になり、小さくするほど幅広く反応するようになります。

大まかな計算式として、検出する周波数の幅（窓の広さ）は以下のように決まります

　検出幅(帯域幅)∽　サンプリング周波数　/　ブロックサイズ

●現在の設定 (7000Hz / 140 = 50Hz): 690Hzを中心に、およそ 665Hz〜715Hz の範囲に強く反応します。
●ブロックサイズを 70 にした場合 (7000Hz / 70 = 100Hz): 700Hzを中心に、およそ 640Hz〜740Hz の幅広い音に反応するようになります（その分、650HzのノイズでもHIGHになりやすくなります）。

速度（時間）と周波数の正確さはトレードオフの関係にあります。

2.3.5 ノイズ対策検

A/D入力の前段にノイズを減らす②バンドパスフィルターは有効です。
デジタル信号処理（ゲーツェル）の前段にアナログフィルターを入れることには、単にノイズを減らすだけでなく、以下の圧倒的な2つのメリットがあります。
a.エイリアシング（折り返し雑音）の防止

現在のサンプリング周波数は 8kHz です。
デジタル信号処理の鉄則として、サンプリング周波数の半分（4kHz）以上の高い周波数成分が入力されると、それが低い周波数のノイズとしてデータ内に「化けて」出現してしまいます（エイリアシング現象）。
前段のフィルターでこれらをカットできます。

b.マイコンのダイナミックレンジの有効活用

今回の無信号時ノイズ（±0.3V）のような不要な成分をはじめに電気回路側で小さく落としておけば、本物の信号に対してADCの分解能（12bit＝4096段階）をフルに贅沢に使うことができます。
これにより、プログラム側の ExclusionSignal（しきい値）を大幅に下げることができ、結果としてより小さな本物の信号でも素早くキャッチできるようになります。

2.4 センター電圧自動取得化

CentralSignalが違っていた（例えば1.533Vを1.2V＝1586とした）場合にはどうなるかと、自動的に中心を出す方法を最終的に説明します。
最初はセンター電圧が変化した場合の影響などから説明します。

2.4.1 センター電圧が変化

　CentralSignal が違っていた場合どうなるか？
　結論から言うと、「非常に低い周波数（直流に近い成分）でブランコが大暴れしてしまい、誤検出や感度低下を引き起こす」ようになります。

(1)なぜ問題になるのか？

本来、CentralSignal は波の中心を 0 にするために引き算しています。
もし、実際の中心が 1.533V (A/D値: 2026) なのに、間違えて 1.2V (A/D値: 1586) を引いてしまうと、引き算が足りず、デジタル空間の波は常に +440 ほどプラス側に浮いた状態になってしまいます。

(2)発生する具体的な不具合

a.無信号時にも power が増大する

ゲーツェルアルゴリズムのループ内で、この浮いた分（+440）が毎回足し算されてしまうため、700Hzの音が鳴っていなくても power の値が非常に大きくなります。
結果として、無信号なのにPB8がHIGHのまま戻らない（誤検出）という現象が起きます。

b.ノイズのしきい値を高くせざるを得ず、感度が落ちる

誤検出を防ぐために ExclusionSignal（しきい値）を高く設定しなければならなくなり、結果として本物の小さな信号を見落とす（感度が悪くなる）ことになります。

2.4.2 センター電圧を出す方法

　自動的に中心を出す方法（自動キャリブレーション）

(1)起動時に「無信号」と仮定して平均を取る（一番簡単）

電源を入れた直後の数十ミリ秒間は信号が入ってこない（静かである）という前提が作れる場合に有効です。
起動時に100回〜200回ほどA/D変換を行い、その平均値を CentralSignal にセットします。

(2)動かしながら「ローパスフィルタ」でジワジワ追従する（高度・確実）

プログラムが動いている最中も、電源電圧の変動などで中心はゆっくりズレることがあります。
そこで、タイマー割り込みの中で「常に最新のA/D値の平均（直流成分）をゆっくり計算し続ける」というデジタルフィルタ（ローパスフィルタ）を入れます。

※この方法を使う場合、CentralSignal の初期値は適当（2048 など）で始めても、電源を入れて数秒で自動的に正しい中心電圧（2026 など）へ勝手に吸い付いて固定されます。

もしシステム的に可能であれば、「方法B（動かしながらジワジワ追従）」を組み込むと、回路の個体差やノイズ環境の変化に非常に強い、頑丈なシステムになります。

2.4.3 センター電圧が変化

前項の(2)で紹介した（動的追従）を組み込んだ完全なプログラム

タイマー割り込み内で CentralSignal をジワジワ更新する処理を追加し、変数の初期化を安全に行うようにした全体のソースコードです。

#include "debug.h"
#include <math.h>

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

// 定数定義
#define SAMPLING_RATE 7000.0f  // サンプリング周波数 (8kHz)
#define BLOCK_SIZE    140      // ブロックサイズ (160サンプル = 20msフレーム)

// ユーザー指定の制御変数
volatile float TargetFreq = 690.0f;          // 600Hz〜800Hzの中心値
volatile int32_t CentralSignal = 2048;      // 【方法B】最初は仮の値(2048)でスタートし、動的に自動追従します
volatile float ExclusionSignal = 110000.0f; // 正規化後のノイズパワー(約78400)より大きく設定

// ゲーツェルアルゴリズム計算用ワーク変数
volatile int sample_count = 0;
float adc_buffer[BLOCK_SIZE];
volatile uint8_t processing_flag = 0;

// GPIO初期化
void GPIO_Config(void) {
    GPIO_InitTypeDef GPIO_InitStructure = {0};
    RCC_APB2PeriphClockCmd(RCC_APB2Periph_GPIOB, ENABLE);

// PB0: アナログ入力 (ADC1_IN8)
    GPIO_InitStructure.GPIO_Pin = GPIO_Pin_0;
    GPIO_InitStructure.GPIO_Mode = GPIO_Mode_AIN;
    GPIO_Init(GPIOB, &GPIO_InitStructure);

// PB8: デジタル出力 (検出フラグ用)
    GPIO_InitStructure.GPIO_Pin = GPIO_Pin_8;
    GPIO_InitStructure.GPIO_Mode = GPIO_Mode_Out_PP;
    GPIO_InitStructure.GPIO_Speed = GPIO_Speed_50MHz;
    GPIO_Init(GPIOB, &GPIO_InitStructure);

GPIO_WriteBit(GPIOB, GPIO_Pin_8, Bit_RESET);
}

// ADC初期化
void ADC_Config(void) {
    ADC_InitTypeDef ADC_InitStructure = {0};
    RCC_APB2PeriphClockCmd(RCC_APB2Periph_ADC1, ENABLE);
    RCC_ADCCLKConfig(RCC_PCLK2_Div6);

ADC_DeInit(ADC1);
    ADC_InitStructure.ADC_Mode = ADC_Mode_Independent;
    ADC_InitStructure.ADC_ScanConvMode = DISABLE;
    ADC_InitStructure.ADC_ContinuousConvMode = DISABLE;
    ADC_InitStructure.ADC_ExternalTrigConv = ADC_ExternalTrigConv_None;
    ADC_InitStructure.ADC_DataAlign = ADC_DataAlign_Right;
    ADC_InitStructure.ADC_NbrOfChannel = 1;
    ADC_Init(ADC1, &ADC_InitStructure);

ADC_Cmd(ADC1, ENABLE);

ADC_ResetCalibration(ADC1);
    while(ADC_GetResetCalibrationStatus(ADC1));
    ADC_StartCalibration(ADC1);
    while(ADC_GetCalibrationStatus(ADC1));
}

// ADC値取得
uint16_t Get_ADC_Val(void) {
    ADC_RegularChannelConfig(ADC1, ADC_Channel_8, 1, ADC_SampleTime_239Cycles5);
    ADC_SoftwareStartConvCmd(ADC1, ENABLE);
    while(!ADC_GetFlagStatus(ADC1, ADC_FLAG_EOC));
    return ADC_GetConversionValue(ADC1);
}

// 8kHzサンプリング用タイマー設定 (TIM2)
void TIM2_Config(void) {
    TIM_TimeBaseInitTypeDef TIM_TimeBaseStructure = {0};
    NVIC_InitTypeDef NVIC_InitStructure = {0};

RCC_APB1PeriphClockCmd(RCC_APB1Periph_TIM2, ENABLE);

// 144MHz動作を想定時の7kHz設定
    TIM_TimeBaseStructure.TIM_Period = 144 - 1;
    TIM_TimeBaseStructure.TIM_Prescaler = 144 - 1;
    TIM_TimeBaseStructure.TIM_ClockDivision = TIM_CKD_DIV1;
    TIM_TimeBaseStructure.TIM_CounterMode = TIM_CounterMode_Up;
    TIM_TimeBaseInit(TIM2, &TIM_TimeBaseStructure);

TIM_ITConfig(TIM2, TIM_IT_Update, ENABLE);

TIM_Cmd(TIM2, ENABLE);
}

// Goertzelフィルター本体の計算関数
uint8_t Process_Goertzel(void) {
    float s_prev = 0.0f;
    float s_prev2 = 0.0f;
    float s;

float omega = (2.0f * M_PI * TargetFreq) / SAMPLING_RATE;
    float coeff = 2.0f * cosf(omega);

// 現在の自動計算されたCentralSignalの値を固定して計算に使用
    float current_center = (float)CentralSignal;

for (int i = 0; i < BLOCK_SIZE; i++) {
        // 自動追従している中心電圧を引き算して直流カット
        float x = adc_buffer[i] - current_center;

s = x + coeff * s_prev - s_prev2;
        s_prev2 = s_prev;
        s_prev = s;
    }

// パワー（振幅の2乗）の計算
    float power = (s_prev * s_prev) + (s_prev2 * s_prev2) - (coeff * s_prev * s_prev2);

// サンプリング数で正規化
    power = power / (BLOCK_SIZE * BLOCK_SIZE);

if (power > ExclusionSignal) {
        return 1; // 検出成功
    } else {
        return 0; // 不検出
    }
}

int main(void) {
    NVIC_PriorityGroupConfig(NVIC_PriorityGroup_2);
    SystemCoreClockUpdate();
    Delay_Init();

GPIO_Config();
    ADC_Config();

// 【初期最適化】起動直後のA/D値を最初の中心電圧としてセットしておく（立ち上がりを早くするため）
    CentralSignal = (int32_t)Get_ADC_Val();

TIM2_Config();

while(1) {
        if (processing_flag) {
            uint8_t detected = Process_Goertzel();

if (detected) {
                GPIO_WriteBit(GPIOB, GPIO_Pin_8, Bit_SET);
            } else {
                GPIO_WriteBit(GPIOB, GPIO_Pin_8, Bit_RESET);
            }
            processing_flag = 0;
        }
    }
}

if (!processing_flag) {
            uint16_t raw_adc = Get_ADC_Val();
            adc_buffer[sample_count] = (float)raw_adc;
            sample_count++;

// 【方法B: 動的ローパスフィルタによる中心電圧の自動追従】
            // 過去の平均値を99.9%残し、新しいデータを0.1%だけブレンドしてジワジワ追いかけます。
            // 計算のオーバーフローや精度低下を防ぐため、1000倍スケールで整数演算しています。
            CentralSignal = (CentralSignal * 999 + raw_adc) / 1000;

if (sample_count >= BLOCK_SIZE) {
                sample_count = 0;
                processing_flag = 1;
            }
        }
    }
}

実装後の動作確認ポイント
電源を入れると、main() の中で最初の1発目のA/D値を読み込んで CentralSignal に代入し、その後はタイマー割り込みの中で自動的に最適な中心（約2026付近）へと勝手に微調整されて落ち着きます。

3.測定　動作(Hi信号出力)値測定

3.1 入力電圧特性

・ExclusionSignal値と入力電圧

限定周波数を検出したことを示すPB8の電圧がHiになる下記可変値を測定する
入力周波数値：690 (Hz)固定
ExclusionSignal値；設定(80,000〜140,000 Step10,000)
入力電圧値:可変(動作時の電圧 RMS)

入力周波数を変化させた時の入力電圧

限定周波数を検出したことを示すPB8の電圧がHiになる下記可変値を測定する
ExclusionSignal値；110,000固定
入力電圧値:設定( 388 , 405 , 421 mV RMS )
入力周波数値：可変(動作時の周波数)

3.2 周波数特性

周波数とプログラム中のpower値の関係を測定する
ExclusionSignal値；110,000固定
入力電圧値:設定 405(mV)RMS固定
入力周波数値：設定(650〜720 Step5Hz )

4.詳細を学ぶための推奨リソース

　ゲーツェルアルゴリズムの数学的な背景や、さらに深い実装ノウハウを知りたい場合は、以下のリソースが非常に参考になります。

4.1 ウェブサイト（日本語で分かりやすいサイト）

4.1.1 DISASSEMBLE CHANNEL（先述のサイト）

ゲーツェルアルゴリズムによる単一周波数検出とDTMF復号
プログラムの実装例や、なぜFFTより軽いのかが視覚的に解説されており、電子工作・マイコン視点で一番わかりやすいです。

4.1.2 Qiita や個人の技術ブログ

検索ワード: ゲーツェルアルゴリズム C言語または Goertzel algorithm dtmf
音声信号処理や通信系のエンジニアが数式付きでコードを解説している記事が多数見つかります。

4.2 書籍（本格的にデジタル信号処理を学びたい場合）

4.2.1 デジタル信号処理の基礎（改訂版）』（著者：岩田彰、コロナ社など）

大学の教科書や技術書で「DFT（離散フーリエ変換）の変形」としてゲーツェルアルゴリズムが紹介されています。

4.2.2 『C言語によるデジタル信号処理入門』（CQ出版・Interface関係の書籍）

マイコンでフィルタや周波数解析を行うための実践的な本です。C言語のコードと数式の結びつきが深く理解できます。

4.3 英語のキーワード（世界水準の情報を探す場合）

より厳密な数学的証明やWikipediaの解説を見る場合は、英語で検索すると良質な情報（シミュレーション用の数式など）が大量に出てきます。

キーワード: Goertzel algorithm Wikipedia / Embedded Goertzel filter tone detection

1.ゲーツェルアルゴリズム

1.1 概略