【トランジスタ技術2006年01月号】　投稿記事『心電計の製作』

Private Const ofsECG As Long = 6 '【Sheet2の格納開始行を示す】

'【変数】
Private TimeCount As Long     '【A/D変換時刻】単位mS
Private aTimeCount As Long    '【現在時刻】
Private bTimeCount As Long    '【過去時刻】
Private aveECG As Long        '【ECG加算平均値】ノッチ処理後
Private servoECG As Long      '【DCサーボ後ECG】
Private CorECG As Long        '【矩形波相関値】
Private RRspan As Long        '【RR間隔】QRS群通しの間隔[mS]
Private HeartRate As Single   '【心拍数】１分間あたりの脈拍数

'【ハムノイズ除去処理関連】
Private pRowInput As Long     '【A/D変換サンプル番号ポインタ】
Private pColTimeCount As Long '【A/D変換時刻管理ポインタ】
Private pColInput As Long     '【A/D変換チャンネルポインタ】
Private Din As Long           '【A/D変換データ】
Private NotchFIFO() As Long   '【NotchFIFOの配列】
Private nNotchFIFO As Long    '【NotchFIFOのサイズ】
Private pNotchFIFO As Long    '【NotchFIFOのポインタ】
Private sNotchFIFO As Long    '【NotchFIFOの累積加算値】
Private kNotchFIFO As Double  '【NotchFIFOの加算平均値を得る為の乗算係数】

'【基線動揺除去処理関連】
Private DCOffset As Long      '【DCオフセット量】

'【矩形波相関処理関連】
Private Const nCorFIFO As Long = 80 '【矩形波相関のサイズ】★必ず4の倍数であること
Private CorFIFO() As Long     '【矩形波相関の配列】
Private pCor1 As Long         '【矩形波相関のポインタ1】初期値　0
Private pCor2 As Long         '【矩形波相関のポインタ2】初期値　nCorFIFO/4
Private pCor3 As Long         '【矩形波相関のポインタ3】初期値　3*nCorFIFO/4
Private sCorFIFO As Long      '【矩形波相関値】

'【QRS群検出時刻処理関連】
Private Const ThreshHold As Long = 40000 '【矩形波相関のピーク検出の為の閾値】
Private indx As Long          '【QRS群検出時刻リストのポインタ】
Private maxCor As Long        '【最大相関値】
Private maxTime As Long       '【最大相関値検出時刻】[mS]
Private PeekList() As Long    '【QRS群検出時刻リストの配列】
Private QRSTimeStamp As Long  '【QRS群検出時刻】

Private Sub CommandButton1_Click()
   pRowInput = 14     '【第14行目から】
   pColTimeCount = 1  '【第1桁目が時刻データ】
   pColInput = 3      '【第2桁目がECGデータ】
   TimeCount = 0      '【時刻番号の初期化】
   aTimeCount = 0     '【A/D変換現在時刻コード】
   bTimeCount = -9999 '【A/D変換前回時刻コード】
   
        InitNotchFIFO    '【初期化】ノッチフィルタ(ハムノイズ除去)
        InitDCServo      '【初期化】DCサーボ処理
        InitRectangleCor '【初期化】矩形波相関
        InitDetectPeek   '【初期化】QRS群検出
LStart:
        Din = GetADvalue()                '【A/D変換値を取得する】
        If (TimeCount < 0) Then End       '【全ての処理を終了】
        aveECG = NotchFilter(Din)         '【ノッチフィルタ処理】
        servoECG = DCServo(aveECG)        '【DCサーボ処理】
        CorECG = RectangleCor(servoECG)   '【矩形波相関処理】
        QRSTimeStamp = DetectPeek(CorECG) '【矩形波相関のピークを検出】
        
        Sheet2.Cells(ofsECG + TimeCount, 1) = TimeCount '【時刻を保存】[mS]
        Sheet2.Cells(ofsECG + TimeCount, 2) = Din       '【ECGのA/D変換値】
        Sheet2.Cells(ofsECG + TimeCount, 3) = aveECG    '【ハムノイズ除去後のECGデータ】
        Sheet2.Cells(ofsECG + TimeCount, 4) = servoECG  '【基線動揺除去後のECGデータ】
        Sheet2.Cells(ofsECG + TimeCount, 5) = CorECG    '【矩形波相関値】
        If QRSTimeStamp > 0 Then '【QRS群が検出された】
           Sheet2.Cells(ofsECG + indx, 6) = QRSTimeStamp - (nCorFIFO \ 2) '【検出時刻】[mS]
        End If
        If indx > 0 Then '【前回検出時刻があれば...】
           RRspan = PeekList(indx) - PeekList(indx - 1)    '【RR間隔】[mS]
           Sheet2.Cells(ofsECG + indx, 7) = RRspan         '
           Sheet2.Cells(ofsECG + indx, 8) = 60000 / RRspan '【心拍数】[回/分]
        End If
        TimeCount = TimeCount + 1 '【時刻更新】
        GoTo LStart
End Sub

Private Sub InitDetectPeek()
   indx = -1            '【QRS群検出時刻を保存するポインタ】
   ReDim PeekList(1000) '【QRS群検出時刻を保存する配列】
   maxCor = 0           '【矩形波相関値の暫定的な最大値】
   maxTime = 0          '【矩形波相関最大値の暫定的な検出時刻】
End Sub

Private Function DetectPeek(Din As Long) As Long
   DetectPeek = -1
   If Din > ThreshHold Then  '【しきい値を超えたか？】
      If maxCor < Din Then   '【前回相関値より大きいか？】
         maxCor = Din        '(Yes) 比較用相関値を保存する
         maxTime = TimeCount '
      Else
         If maxTime + 1 = TimeCount Then
            DetectPeek = TimeCount   '【関数出力】
            indx = indx + 1          '【リスト保存ポインタを更新する】
            PeekList(indx) = maxTime 'ピークを越えた⇒時刻をリストに保存
         End If
      End If
   Else
      maxCor = 0 '【閾値以下】
   End If
End Function

Private Sub InitRectangleCor()
   Dim i As Long
   
   ReDim CorFIFO(nCorFIFO - 1) '【矩形波相関処理用FIFO配列初期化】
   For i = 0 To nCorFIFO - 1   '
      CorFIFO(i) = 0           '
   Next i                      '
   pCor1 = 3 * nCorFIFO \ 4    '【管理ポインタの初期化】
   pCor2 = nCorFIFO \ 4        '
   pCor3 = 0                   '
   sCorFIFO = 0                '【矩形波相関値】
End Sub

Private Function RectangleCor(Din As Long) As Long
   sCorFIFO = sCorFIFO - Din + 2 * CorFIFO(pCor1) - 2 * CorFIFO(pCor2) + CorFIFO(pCor3) '【矩形波相関演算】
   CorFIFO(pCor3) = Din                                                                 '【新規データをFIFOに保存】
   pCor1 = pCor1 + 1                  '【ポインタpCor1の更新】
   If pCor1 = nCorFIFO Then pCor1 = 0 '
   pCor2 = pCor2 + 1                  '【ポインタpCor2の更新】
   If pCor2 = nCorFIFO Then pCor2 = 0 '
   pCor3 = pCor3 + 1                  '【ポインタpCor3の更新】
   If pCor3 = nCorFIFO Then pCor3 = 0 '
   RectangleCor = sCorFIFO            '【矩形波相関値を出力する】
End Function

Private Sub InitDCServo()
   DCOffset = 0 '【DCオフセットをリセットする】
End Sub

Private Function DCServo(Din As Long) As Double
   DCOffset = DCOffset + (Din - DCOffset) * 0.01 '【現DCオフセットとの差分値を累積加算してDCオフセットを得る】
   DCServo = Din - DCOffset                      '【DCオフセットを差し引いた分を出力する】
End Function

Private Sub InitNotchFIFO() '【ハムノイズ除去用FIFOのサイズを決め、0クリアする】
   Dim i As Long
   
   If (OptionButton1.Value = True) Then nNotchFIFO = 20 Else nNotchFIFO = 17
   ReDim NotchFIFO(nNotchFIFO - 1) '【配列NotchFIFOを確保】
   For i = 0 To nNotchFIFO - 1     '【配列NotchFIFOを0クリア】
      NotchFIFO(i) = 0             '
   Next i                          '
   pNotchFIFO = 0                  '【ポインタの初期化】
   sNotchFIFO = 0#                 '【NotchFIFOの累積加算値】
   kNotchFIFO = 1# / nNotchFIFO    '【加算平均値を得るための乗算係数】
End Sub

Private Function NotchFilter(Din As Long) As Double
   sNotchFIFO = sNotchFIFO + Din - NotchFIFO(pNotchFIFO) '【累積加算値を求める】
   NotchFIFO(pNotchFIFO) = Din                           '【FIFOに新規データを格納】
   pNotchFIFO = pNotchFIFO + 1                           '【ポインタ修正】
   If (pNotchFIFO = nNotchFIFO) Then pNotchFIFO = 0      '【ポインタ修正】
   NotchFilter = sNotchFIFO * kNotchFIFO                 '【加算平均値を返す】
End Function

Private Function GetADvalue() As Long
   aTimeCount = Sheet1.Cells(pRowInput, pColTimeCount) '【現在時刻情報を読込む】
   If (aTimeCount < bTimeCount) Then                   '【過去値と比較】⇒減少なら測定終了
      TimeCount = -1 '【データ終了】
      GetADvalue = 0
   Else
      GetADvalue = Sheet1.Cells(pRowInput, pColInput)  '【A/D変換値を読込む】
      pRowInput = pRowInput + 1                        '【読込ポインタを進める】
      bTimeCount = aTimeCount                          '【現在値⇒過去値】
   End If
End Function