【光学系】波長から色座標を求める

やること

　適当な測定波長などから、CIE1931色空間のxy色度を求めてみる。
　今回、すべてCopilotに任せてみた。

結果

　以下の通り。

CreateColorMatchingTable_AndChart

Sub CreateColorMatchingTable_AndChart()

    Dim ws As Worksheet
    Set ws = ActiveSheet

    Dim λ As Long
    Dim row As Long
    row = 2

    '--- 見出し ---
    ws.Range("A1:H1").Value = Array("Wavelength (nm)", "X", "Y", "Z", _
                                    "Measured", "x̅(λ)", "y̅(λ)", "z̅(λ)")

    '--- 波長 400〜700 nm を 10nm 刻みで入力 ---
    For λ = 400 To 700 Step 10
        ws.Cells(row, 1).Value = λ
        row = row + 1
    Next λ

    '--- 等色関数データ（400〜700nm） ---
    Dim cmfX, cmfY, cmfZ

    cmfX = Array(0.001, 0.006, 0.023, 0.059, 0.087, 0.096, 0.088, 0.063, 0.032, 0.011, _
                 0.002, 0.003, 0.020, 0.052, 0.094, 0.148, 0.208, 0.261, 0.315, 0.390, _
                 0.447, 0.449, 0.399, 0.308, 0.219, 0.143, 0.084, 0.044, 0.023, 0.011, 0.005)

    cmfY = Array(0.000, 0.000, 0.001, 0.002, 0.006, 0.011, 0.018, 0.029, 0.047, 0.072, _
                 0.109, 0.163, 0.221, 0.269, 0.309, 0.339, 0.349, 0.326, 0.299, 0.288, _
                 0.265, 0.225, 0.178, 0.127, 0.086, 0.054, 0.031, 0.016, 0.008, 0.004, 0.002)

    cmfZ = Array(0.007, 0.028, 0.109, 0.287, 0.439, 0.504, 0.508, 0.416, 0.275, 0.160, _
                 0.092, 0.051, 0.024, 0.013, 0.007, 0.003, 0.001, 0.001, 0.001, 0.000, _
                 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000)

    '--- 等色関数を F〜H 列に書き込む ---
    row = 2
    Dim i As Long
    For i = 0 To UBound(cmfX)
        ws.Cells(row, 6).Value = cmfX(i)
        ws.Cells(row, 7).Value = cmfY(i)
        ws.Cells(row, 8).Value = cmfZ(i)
        row = row + 1
    Next i

    '--- 白色LEDの測定値（簡易モデル） ---
    row = 2
    For i = 400 To 700 Step 10
        Dim v As Double
        v = Exp(-((i - 450) ^ 2) / (2 * 12 ^ 2)) * 1.0
        v = v + Exp(-((i - 560) ^ 2) / (2 * 40 ^ 2)) * 0.6
        ws.Cells(row, 5).Value = v
        row = row + 1
    Next i

    '--- B〜D列に X,Y,Z の計算式を入れる ---
    row = 2
    Do While ws.Cells(row, 1).Value <> ""
        ws.Cells(row, 2).Formula = "=E" & row & "*F" & row
        ws.Cells(row, 3).Formula = "=E" & row & "*G" & row
        ws.Cells(row, 4).Formula = "=E" & row & "*H" & row
        row = row + 1
    Loop

    '===============================
    '   グラフ作成（波長を横軸）
    '===============================
    Dim cht As ChartObject
    Dim lastRow As Long
    lastRow = ws.Cells(ws.Rows.Count, "A").End(xlUp).Row

    For Each cht In ws.ChartObjects
        cht.Delete
    Next cht

    Set cht = ws.ChartObjects.Add(Left:=350, Top:=20, Width:=650, Height:=380)
    cht.Chart.ChartType = xlLine

    With cht.Chart
        '--- 等色関数（点線） ---
        .SeriesCollection.NewSeries
        .SeriesCollection(1).Name = "x̅(λ)"
        .SeriesCollection(1).Values = ws.Range("F2:F" & lastRow)
        .SeriesCollection(1).XValues = ws.Range("A2:A" & lastRow)
        .SeriesCollection(1).Format.Line.ForeColor.RGB = RGB(255, 0, 0)
        .SeriesCollection(1).Format.Line.DashStyle = msoLineDash

        .SeriesCollection.NewSeries
        .SeriesCollection(2).Name = "y̅(λ)"
        .SeriesCollection(2).Values = ws.Range("G2:G" & lastRow)
        .SeriesCollection(2).XValues = ws.Range("A2:A" & lastRow)
        .SeriesCollection(2).Format.Line.ForeColor.RGB = RGB(0, 180, 0)
        .SeriesCollection(2).Format.Line.DashStyle = msoLineDash

        .SeriesCollection.NewSeries
        .SeriesCollection(3).Name = "z̅(λ)"
        .SeriesCollection(3).Values = ws.Range("H2:H" & lastRow)
        .SeriesCollection(3).XValues = ws.Range("A2:A" & lastRow)
        .SeriesCollection(3).Format.Line.ForeColor.RGB = RGB(0, 0, 255)
        .SeriesCollection(3).Format.Line.DashStyle = msoLineDash

        '--- XYZ（実測×等色関数） ---
        .SeriesCollection.NewSeries
        .SeriesCollection(4).Name = "X"
        .SeriesCollection(4).Values = ws.Range("B2:B" & lastRow)
        .SeriesCollection(4).XValues = ws.Range("A2:A" & lastRow)
        .SeriesCollection(4).Format.Line.ForeColor.RGB = RGB(255, 0, 0)

        .SeriesCollection.NewSeries
        .SeriesCollection(5).Name = "Y"
        .SeriesCollection(5).Values = ws.Range("C2:C" & lastRow)
        .SeriesCollection(5).XValues = ws.Range("A2:A" & lastRow)
        .SeriesCollection(5).Format.Line.ForeColor.RGB = RGB(0, 180, 0)

        .SeriesCollection.NewSeries
        .SeriesCollection(6).Name = "Z"
        .SeriesCollection(6).Values = ws.Range("D2:D" & lastRow)
        .SeriesCollection(6).XValues = ws.Range("A2:A" & lastRow)
        .SeriesCollection(6).Format.Line.ForeColor.RGB = RGB(0, 0, 255)

        '--- 測定値（黒） ---
        .SeriesCollection.NewSeries
        .SeriesCollection(7).Name = "Measured"
        .SeriesCollection(7).Values = ws.Range("E2:E" & lastRow)
        .SeriesCollection(7).XValues = ws.Range("A2:A" & lastRow)
        .SeriesCollection(7).Format.Line.ForeColor.RGB = RGB(0, 0, 0)

        .HasTitle = True
        .ChartTitle.Text = "Color Matching Functions & White LED Spectrum"
        .Axes(xlCategory).HasTitle = True
        .Axes(xlCategory).AxisTitle.Text = "Wavelength (nm)"
        .Axes(xlValue).HasTitle = True
        .Axes(xlValue).AxisTitle.Text = "Value"
    End With

    '===============================
    '   合計（B33〜D33）と xy（B34, C34）をグラフ作成後に配置
    '===============================
    ws.Range("A34").Value = "(x,y)="

    ws.Range("B33").Formula = "=SUM(B2:B32)"
    ws.Range("C33").Formula = "=SUM(C2:C32)"
    ws.Range("D33").Formula = "=SUM(D2:D32)"

    ws.Range("B34").Formula = "=B33/(B33+C33+D33)"
    ws.Range("C34").Formula = "=C33/(B33+C33+D33)"

    '--- 罫線を追加（A33〜D34） ---
    With ws.Range("A33:D34").Borders
        .LineStyle = xlContinuous
        .Weight = xlThin
    End With

    With ws.Range("A33:D34").Borders(xlEdgeLeft)
        .LineStyle = xlContinuous
        .Weight = xlMedium
    End With
    With ws.Range("A33:D34").Borders(xlEdgeTop)
        .LineStyle = xlContinuous
        .Weight = xlMedium
    End With
    With ws.Range("A33:D34").Borders(xlEdgeRight)
        .LineStyle = xlContinuous
        .Weight = xlMedium
    End With
    With ws.Range("A33:D34").Borders(xlEdgeBottom)
        .LineStyle = xlContinuous
        .Weight = xlMedium
    End With

    MsgBox "表・合計・xy・グラフの作成が完了しました！"

End Sub

補足

等色関数データ（400〜700nm）

色を科学する　＜番外編＞ XYZを自分で計算してみよう！｜［色彩工学/カラーサイエンス］YoshiColor
　ここの等色関数の表のデータを転記したようです。ゆえに400〜700nmの10nm刻みの簡略版であり、CIEの公式データをそのまま使ったものではありません。実際にはJIS辺りを参照して頂くのがよいかと思います。

日本産業標準調査会：データベース検索-JIS検索

(追記) JIS Z8781-1(2012)

　ここの等色関数の表のデータを転記したようです。ゆえに400〜700nmの10nm刻みの簡略版であり、CIEの公式データをそのまま使ったものではありません。実際にはJIS辺りを参照して頂くのがよいかと思います。

　「表1-CIE1931 測色標準観測者の等色関数及び色度座標」の波長400nmの値と比較。

今回	JIS
0.001	0.014 310 00
0.000	0.000 396 000 0
0.007	0.067 850 01

　同様に波長460nmの値と比較。

今回	JIS
0.088	0.290 800 0
0.018	0.060 000 00
0.508	1.669 200 0

三刺激値は，一般的に相対値で評価されるので，

と、JISに記載がありますので、この辺りはJISを各自ご確認ください。

白色LEDの測定値（簡易モデル）

これは 実測データの引用ではなく、白色LEDの典型的な形状を模倣した合成スペクトル(2つのガウス分布を足し合わせたもの)です。
したがって、特定のメーカーや論文のデータではありません。

青ピーク（450nm付近）

白色LEDの励起光（青色LED）を模倣

中心波長：450 nm
半値幅：狭い（σ=12）
強度：1.0

黄色蛍光体の広帯域発光（560nm付近）

白色LEDの蛍光体（YAG系など）を模倣

中心波長：560 nm
半値幅：広い（σ=40）
強度：0.6

まとめ

色空間とは？開発に役立つ使い方、トレンド記事やtips - Qiita
　正確な話は、ここ辺りを参照頂きたい。個人的に以下をメモしておく。

CIE 関連の色空間について #色空間 - Qiita
投稿日 2025年07月02日