\documentclass[dvipdfmx]{jsarticle}
\usepackage{graphicx}
\usepackage{listings,jlisting}
\usepackage{ice-report-cover}
\usepackage{fancyhdr}
\usepackage{lastpage}
\usepackage{physics}
\usepackage{url}
\usepackage{mathtools}
\usepackage{amsmath,amssymb}
\usepackage[binary-units=true]{siunitx}
\AtBeginDocument{\RenewCommandCopy\qty\SI}
\usepackage{float}
\usepackage{tikz}
\usepackage{circuitikz}
\usepackage[hang,small,bf]{caption}
\usepackage[subrefformat=parens]{subcaption}

\lstset{
  basicstyle={\ttfamily},
  identifierstyle={\small},
  commentstyle={\smallitshape},
  keywordstyle={\small\bfseries},
  ndkeywordstyle={\small},
  stringstyle={\small\ttfamily},
  captionpos=b,
  frame={tb},
  breaklines=true,
  columns=[l]{fullflexible},
  numbers=left,
  xrightmargin=0zw,
  xleftmargin=3zw,
  numberstyle={\scriptsize},
  stepnumber=1,
  numbersep=1zw,
  lineskip=-0.5ex
}
\renewcommand{\lstlistingname}{リスト}

\fancypagestyle{mypagestyle}{
  \lhead{}
  \chead{}
  \rhead{\thepage/\pageref{LastPage}}%ヘッダ左を空に
  \cfoot{}
  \renewcommand{\headrulewidth}{0.0pt}%ヘッダの線を消す
}

\renewcommand{\refname}{}
\makeatletter
\renewcommand{\@cite}[1]{\raisebox{1ex}{\footnotesize \,#1)}}
\renewcommand{\@biblabel}[1]{#1)}
\makeatother

\報告者{5}{}{G}{22330}{野口 煌陽}
\共同実験者*{}{}{}{}{}{}{}
\プレ受付無し

\begin{document}
\実験課題{4}{音声データの信号処理}
\実験場所{情報処理}
\実験日{2026/4/20}{2026/4/27}
\実験条件{}{}{}{}
\担当教職員{川本}
\提出期限{2026/5/11}
\提出日{\the\year/\the\month/\the\day}
\coverprint
\selfchecksheet

\newpage
\pagestyle{empty}
\begin{abstract}
本実験は音声データの分析を対象に，1次元の時系列信号に対する処理の実装を通じて，信号処理の基礎を理解することを目的としている．
離散フーリエ変換(DFT)とは，標本化されたディジタル信号に対する周波数解析を行うための方法である．
本実験では，DFTおよびDFTを高速化した高速フーリエ変換(FFT)をC言語で実装し，それらの実行効率や計算結果の比較を行った．
さらに，窓関数を適用するプログラムを実装し，振幅スペクトルにどのような変化が現れるかについての確認と考察を行った．
分析する音声データには，標本化周波数\qty{16}{k\hertz}，量子化ビット数\qty{16}{\bit}，チャネル数1chのものを使用した． 
各実験課題の実験結果をgnuplotを用いてグラフとして描画し，それをもとに実験結果を考察した．
振幅スペクトルの計算と可視化の実験では，どちらのフーリエ変換で周波数解析を行っても結果は変わらないが，DFTの方がFFTに比べ約345倍もの実行時間の差があることが確認された．
窓関数の導入の実験では，窓関数により振幅スペクトルが全体的に小さくなることが確認された．
各母音の音声の周波数解析の実験では，発声内容によって振幅スペクトルのピーク値が異なることが確認された．
さらに発展課題として，各自の連続した音声データ全体のスペクトログラムを生成し，各母音のフォルマントの違いや，連続発声時における声道形状の変化に伴う動的な周波数遷移の様子を可視化・分析した．
これらの実験を通じて，離散フーリエ変換による音声データの周波数解析の動作原理と実装方法を総合的に理解することができた．
また，考察事項として窓関数適用による振幅値の減少，ダウンサンプリングの副作用，窓関数と振幅スペクトルの関係についての考察を行った．
以下にその要点を述べる：
考察1）窓関数の値域が$0\leqq w[n]\leqq1$であるため，振幅スペクトルは小さくなる．本実験では，信号の平均パワーが窓関数適用により約0.27倍になることが確認された．
考察2）前処理をせず，単純に観測信号を間引くだけでは高音域の情報の損失や折り返し雑音が発生してしまう．折り返し雑音対策のために，前処理としてアンチエイリアシングフィルタを適用する．
考察3）窓関数を特徴づける要素として，メインローブの幅とサイドローブのレベルがある．これらの違いにより，様々な窓関数が存在するため，状況に応じて選択する必要がある．

\end{abstract}
\newpage
\setcounter{page}{1}
\pagestyle{mypagestyle}
\section{目的}
音声データの分析を対象に，1次元の時系列信号に対する処理の実装を通じて，信号処理の基礎を理解する．
\section{原理}
\subsection{アナログ信号とディジタル信号}
音声は，空気の振動などにより伝播する．
音声を信号として捉えると，振幅が時間的に変化することで情報が伝わる．
音声の情報は，マイクロフォンにより電気信号として観測できる．
このように観測される信号は，時間変化する1次元の情報であり，音声波形として時間成分と振幅成分として表現できる．
電気信号の時点では，時間成分及び振幅成分ともに連続値であり，「アナログ信号」という．

音声信号をコンピュータに取り込み，処理するためには，時間成分及び振幅成分ともに離散値に変換した「ディジタル信号」にする必要がある．
アナログ信号をディジタル信号に変換することをA/D変換という．
逆に，ディジタル信号をアナログ信号に変換することをD/A変換という．
コンピュータ上の音声信号をスピーカーやヘッドフォンで再生する際にはD/A変換が必要となる．

\begin{figure}[h]
  \centering
  \includegraphics[scale=0.25]{ADconvFlow.png}
  \caption{A/D変換の流れ} \label{fig:ADconvFlow}
\end{figure}

入力のアナログ信号をコンピュータに取り込むディジタル信号に変換するための基本的な処理の流れを図\ref{fig:ADconvFlow}に示す．
A/D変換は2種の離散化により実現される．
時間方向に離散化することを標本化（サンプリング）といい，振幅方向に離散化することを量子化という．

標本化は，アナログ信号から一定の時間間隔毎に信号値を得る操作である．
信号値を得る時間間隔$T$を標本化間隔あるいは標本化周期，1秒あたりの標本化回数($=\frac{1}{T}$)を標本化周波数という．
標本化定理\footnote{信号の最大周波数が$f_{max}[\mathrm{Hz}]$のとき，信号を$2f_{max}[\mathrm{Hz}]$より高い標本化周波数で標本化すると元の信号を復元することができる．}により，標本化対象の信号は，標本化周波数の半分未満に帯域制限されている必要がある．
このため，標本化前にはアンチエイリアスフィルタと呼ばれる低域通過フィルタを通し，不要な高域成分を除去する．
もし，信号の最大周波数の2倍未満である低い周波数で標本化を行なうと，見かけ上，低い周波数が存在しているかのように見える現象が生じる．このような現象をエイリアシング（折り返し歪）という．

量子化は，ある範囲の振幅値を代表値に置き換えることで離散的な値に変換する．
このため，量子化では代表値に置き換える際に誤差を生じることとなる．
この誤差を量子化誤差という．
また，量子化した振幅値を何ビットの整数値で表現するかにより，振幅情報をなめらかさが変化する．この情報を量子化ビット数という．
量子化ビット数が大きいほど，多くの代表値で表現することでなめらかな振幅が表現でき，量子化誤差が小さくなる．
本実験では，波形符号化方式として線形パルス符号化(linear pulse code modulation: linear PCM)を採用する．
これは，アナログ信号を一様なステップで量子化するもので，通常A/D変換として行われているものと同じである．

\subsection{離散フーリエ変換}
標本化定理を満たすように標本化周期$T$で標本化されたディジタル信号$x[n]$について考える．
ディジタル信号に対する周波数解析を行うには，離散化されたフーリエ変換，すなわち離散フーリエ変換(DFT;discrete Fourier transform)が用いられる．

周期$N$点の周期信号と見なし，$N$点のディジタル信号に対する離散フーリエ変換は式\eqref{eq:dft}のように定義される．
\begin{equation}
  \label{eq:dft}
  X[k]=\sum_{n=0}^{N-1}x[n]e^{-j\frac{2\pi}{N}kn}=\sum_{n=0}^{N-1}x[n]W_N^{kn}
\end{equation}
また，逆離散フーリエ変換は式\eqref{eq:idft}のように定義される．
\begin{equation}
  \label{eq:idft}
  x[n]=\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j\frac{2\pi}{N}kn}=\sum_{k=0}^{N-1}X[k]W_N^{-kn}
\end{equation}
ただし，$W_N=e^{-j\frac{2\pi}{N}}$，$n$は時間に関する整数変数，$k$は周波数に関する整数変数である．
$n$は切り出した$N$点のディジタル信号の先頭を0としたときの時間$t=nT$と対応する．
$k$は各周波数$\omega=\frac{2\pi}{T}\frac{k}{N}$に対応する．
なお，$N$点のDFTの時間計算量は$O(N^2)$となる．

\subsection{高速フーリエ変換}
DFTでは分析点数$N$の2乗オーダーの演算量を必要とする．
このため，分析点数が多い場合や繰り返しフーリエ変換を行う場合には，その演算量が問題となることも多い．
そのため，多くの場合で高速フーリエ変換(FFT;fast Fourier transform)という方法が用いられる．
ここでは分析点数$N$が2のべき乗の形で与えられる条件下で，時間分割法に基づくFFTを紹介する．

時間分割法では，分析対象の信号を「奇数番目のデータ」と「偶数番目のデータ」に分けていくことで，問題を小分けし，計算量を削減する．
ここでは分析点数$N=8$を例として考える．

まずは信号$x[n]$を，
\begin{align*}
\begin{cases}
  偶数番目の信号\quad x_0[n]=x[2n] & (n=0,1,2,3) \\
  奇数番目の信号\quad x_1[n]=x[2n+1]\quad & (n=0,1,2,3)
\end{cases}
\end{align*}
に分割する．
このとき，信号$x[n]$のDFTは，$W_N=e^{-j\frac{2\pi}{N}}$とおくと，式\eqref{eq:dft_divide}のように表現できる．
\begin{align}
  X[k]&=\sum_{n=0}^7x[n]W_8^{kn} \notag\\
      &=\sum_{n:even}x[n]W_8^{kn}+\sum_{n:odd}x[n]W_8^{kn} \notag\\
      &=\sum_{n=0}^3x_0[n]W_8^{k(2n)}+\sum_{n=0}^3x_1[n]W_8^{k(2n+1)} \notag\\
      &=\underbrace{\sum_{n=0}^3x_0[n]W_4^{kn}}_{x_0[n]\text{のDFT}}+W_8^k\underbrace{\sum_{n=0}^3x_1[n]W_4^{kn}}_{x_1[n]\text{のDFT}} \label{eq:dft_divide}
\end{align}
つまり，$x_0[n]$のDFTを$X_0[k]$，$x_1[n]$のDFTを$X_1[k]$とおくと，
\[
  X[k]=X_0[k]+W_8^kX_1[k]\quad(k=0,1,2,3)
\]
となる．
ここで$k$は0,1,2,3であり，8点のDFTの前半4点分の計算のみである．
同様に後半4点分の計算$X[k+4]$について考える．
式\eqref{eq:dft_divide}において，$k$に$k+4$を代入すると，次のようになる．
\begin{equation}
  \label{eq:fft_k+4}
  X[k+4]=\sum_{n=0}^3x_0[n]W_4^{(k+4)n}+W_8^{k+4}\sum_{n=0}^3x_1[n]W_4^{(k+4)n}
\end{equation}
ここで，$W_N=e^{-j\frac{2\pi}{N}}$を考慮すると，
\begin{align}
  W_4^{(k+4)n}&=W_4^{kn} \\
  W_8^{(k+4)}&=-W_8^k
\end{align}
となるため，これらを式\eqref{eq:fft_k+4}に代入すると，
\begin{align}
  X[k+4]&=\sum_{n=0}^3x_0[n]W_4^{kn}-W_8^k\sum_{n=0}^3x_1[n]W_4^{kn} \notag \\
        &=X_0[k]-W_8^kX_1[k]
\end{align}
となる．
以上より，
\begin{align*}
\begin{cases}
  X[k] &= X_0[k]+W_8^kX_1[k]\quad (n=0,1,2,3) \\
  X[k+4] &= X_0[k]-W_8^kX_1[k]\quad (n=0,1,2,3) 
\end{cases}  
\end{align*}
を得る．

4点のDFTについても同様に2点のDFTに分割することができる．
このように最終的に2点のDFTまで分割し，計算結果を統合することでFFTが実現される．
なお，$N$点間の時間分割法によるFFTの時間計算量は$O(N\log N)$となる．

\subsection{信号の周波数分析}
信号$x[n]$の離散フーリエ変換で得られるスペクトル$X[k]$は複素数であり，図\ref{fig:spectrum}のように実数成分$Re(X[k])$と虚数成分$Im(X[k])$から，振幅成分である振幅スペクトル$\abs{X[k]}$と，位相成分である位相スペクトル$\angle X[k]$に分解できる．
\begin{align*}
\begin{cases}
  \ 振幅スペクトル\quad \abs{X[k]}\quad&=\quad\sqrt{\{Re(X[k])\}^2+\{Im(X[k])\}^2}\\
  \ 位相スペクトル\quad \angle X[k]\quad&=\quad\tan^{-1}\left(\frac{Im(X[k])}{Re(X[k])}\right)
\end{cases}
\end{align*}
\begin{figure}[h]
  \centering
  \includegraphics[scale=0.25]{spectrum.png}
  \caption{振幅スペクトルと位相スペクトル} \label{fig:spectrum}
\end{figure}

\subsection{窓関数}
DFTやFFTでは，分析点数$N$を周期とする信号と見なして分析される．
そのため，切り出し区間の両端点で異なる値を持つ場合，信号における不連続点が生じてしまい，元の信号に含まれない周波数成分が現れることがある．
この現象をスペクトル漏れという．
このため，切り出し区間の両端点付近を0に近くなるような関数を掛け合わせることで，両端点での波形の急激な変化を抑える方法が用いられる．
このように，信号を指定区間に切り出すための関数を窓関数という．
元の信号を$x[n]$，窓関数を$w[n]$とすると，窓関数により切り出された信号$y[n]$は次のように表される．
\[
  y[n]=x[n]\cdot w[n]
\]

窓関数としては，用途に応じて様々なものが存在する．
例えば，図\ref{fig:hanning}に示すハニング窓は，以下の式で表現される．
\begin{equation}
  w[n]=0.5-0.5\cos\left(\frac{2\pi n}{N-1}\right)
  \label{eq:hanning}
\end{equation}
\begin{figure}[h]
  \centering
  \includegraphics[scale=0.4]{hanning.pdf}
  \caption{ハニング窓の例} \label{fig:hanning}
\end{figure}

\subsection{音声信号の周波数分析}
音声信号は，信号に含まれる周波数特性を時間変化させることにより情報を伝達している．
このため，音声の周波数分析を行う場合は，時間変化する信号であっても「切り出した短時間の中では周波数特性が変化しない」と仮定し，切り出した短時間の信号を周期信号と見なしてDFTやFFTを適用することで，周波数分析を行う．
このような短時間の信号に対する分析を，一定周期で切り出す区間をずらしながら行う．
このとき，短時間の信号に切り出す処理をフレーム化処理といい，切り出した信号をフレーム，フレームに含まれる信号点数をフレーム長，フレームをずらす時間間隔をフレーム周期という，
フレーム化処理のイメージ図を図\ref{fig:frame}に示す．
\begin{figure}[h]
  \centering
  \includegraphics[scale=0.25]{frame.png}
  \caption{フレーム化処理のイメージ図} \label{fig:frame}
\end{figure}
フレーム化処理により得られた各フレームに対し，必要に応じて窓関数をかけ，FFTを適用した結果から振幅スペクトルを求めることで，
フレームを切り出した時刻における振幅スペクトルを推定することができる．

フレーム毎に得られる振幅スペクトルを時系列と見なし，時間・周波数・信号成分の3次元の情報として表現したものをスペクトログラムという．
一般的なスペクトログラムの可視化方法としては，時間を横軸，周波数を縦軸，振幅を輝度や色で表現した2次元画像を用いることが多い．

振幅スペクトル$\abs{X[k]}$やパワースペクトル$\abs{X[k]}^2$を扱うとき，対数で表現することも多い．
音声においても信号の比として相対量を可視化することが多い．
この時用いられるのがデシベル(dB)単位での表現であり，振幅スペクトルの場合は，以下のように表される．
\[
20\log_{10}\frac{\abs{X[k]}}{\abs{X_b}}\ [\mathrm{dB}]
\]

ここで$X_b$は振幅スペクトルの相対量を表現するための基準値であり，可視化などにおいては「最大振幅値」を基準値とすることが多い．
つまり，最大振幅値を0$[\mathrm{dB}]$とし，可視化のために設定する下限値\footnote{対数において$\log 0\to -\infty$となることもあり，可視化の目的に沿ってユーザが設定したある一定の幅で可視化する．}までの幅でスペクトルあるいはスペクトログラムを可視化する．
例えば信号の振幅値の量子化ビット数が16$\mathrm{bit}$(65536段階)である場合，信号のダイナミックレンジ\footnote{ダイナミックレンジは，識別可能な信号の最小値と最大値の比を表す情報のことである．}は，
\[
20\log_{10}65536 \simeq 96[\mathrm{dB}]
\]
となる．
なお符号付き2バイト整数で振幅値を扱う場合，最大振幅値は32767であり，実質的な信号のダイナミックレンジ90$[\mathrm{dB}]$程度となる．
そのため，振幅スペクトルの最大値を0$[\mathrm{dB}]$としたとき，$0\sim -90[\mathrm{dB}]$程度の範囲で十分な情報を可視化できる．
実際には，観測した信号の最大振幅などによって調整する必要がある．

\section{実験課題}
以下の課題を実現するC言語のプログラム，およびgnuplotスクリプトを作成するとともに，所望の図を得るまでの手順を示せ．
なお，信号処理のための各プログラムはC言語を用いて実装し，C言語の仕様に付随して提供される標準的な関数のみを使用して作成すること．
また，原理のアルゴリズムを対応付けてプログラムの概要を説明すること．

\subsection*{課題1 バイナリファイルの読み込み}
ヘッダのない音声データsample01.swは以下の条件でディジタル化された音声波形のデータである：
\begin{description}
  \item 標本化周波数：\qty{16}{k\hertz}
  \item 量子化ビット数：\qty{16}{\bit}
  \item チャネル数：\qty{1}{ch}(モノラル)
  \item 符号化方式：Linear PCM
\end{description}
1サンプル分のデータは\qty{2}{bytes}で表現され，C言語における，signed short型の整数変数1つに対応する．
この信号の先頭30ミリ秒(msec)分のデータを読み込み，波形データとしてgnuplotを用いて可視化せよ．
なお，縦軸は振幅値（単位なし），横軸は時間$[\mathrm{msec}]$として可視化すること．
また，読み込む際，signed short型変数で期待される最大値において，最大振幅の絶対値が1となるように振幅を正規化した実数値として扱うこと（つまり，サンプル値32767が振幅値1となるように正規化）．

\subsection*{課題2 $DFT$および$FFT$による振幅スペクトルの計算と可視化}
音声データsample01.swについて，先頭512点を読み込み，フーリエ変換することで振幅スペクトルを算出し，$0\sim 8$\si{k\hertz}の範囲で振幅スペクトルを可視化せよ．
なお，縦軸は振幅スペクトル値（単位なし），横軸は周波数[$\mathrm{Hz}$]として可視化すること．
また，同じデータについて縦軸振幅スペクトル値[$\mathrm{dB}$]，横軸は周波数[$\mathrm{Hz}$]として可視化する場合についても合わせて示すこと．
フーリエ変換については離散フーリエ変換（DFT;Discrete Fourier Transform）を実現する関数dftおよび高速フーリエ変換（FFT;Fast Fourier Transform）を実現する関数fftを実装し，その動作原理と使い方，計算結果の比較，計算速度の比較をレポートに示すこと．

\subsection*{課題3 窓関数の導入}
音声データsample01.swについて，先頭512点をそのまま読み込んだ後，ハニング窓（hanning window）を適用する．
窓長$N$は512となる．
観測信号$x[n]$と窓関数適用後の信号$y[n]$をそれぞれフーリエ変換し，振幅スペクトルを可視化することで比較せよ．
なお，縦軸は振幅スペクトル（単位なし），横軸は周波数[$\mathrm{Hz}$]とし，横軸の範囲$0\sim 8$\si{k\hertz}にて可視化すること．
また，同じデータについて縦軸は振幅スペクトル[$\mathrm{dB}$]，横軸は周波数[$\mathrm{Hz}$]として可視化する場合についても合わせて示すこと．

\subsection*{課題4 各自の音声の周波数分析}
レポート作成者本人の音声を対象として振幅スペクトルを計算・可視化せよ．
なお，適切な部分を複数個所切り出し，それぞれの振幅スペクトルを比較することで，発生内容による振幅スペクトルの違いを分析すること．

\subsection*{発展課題1 各自の音声のスペクトログラムの可視化}
レポート作成者本人の音声を対象として，周波数分析を行い，スペクトログラムを可視化せよ．
スペクトログラムはモノクロ表示とし，振幅スペクトル値[$\mathrm{dB}$]が大きいところを黒，小さいところを白とするように可視化すること．
可視化の形式としては，ヒートマップを参考にするとよい．
スペクトログラムは縦軸を周波数(0～8\si{k\hertz})，横軸を時間[$\mathrm{sec}$]とすること．
また，余裕があればモノクロのスペクトログラムと合わせて，カラーのスペクトログラムの可視化も示すと良い．
必ず，複数の音声について分析結果を掲載し，その違いを分析すること．

なお，音声の分析条件は以下に従うものとする．
\begin{description}
  \item 標本化周波数：\qty{16}{k\hertz}
  \item 量子化ビット数：\qty{16}{\bit}（符号化方式：Linear PCM）
  \item チャネル数： \qty{1}{ch}（モノラル）
  \item フレーム長（窓長）：\qty{32}{m\sec}（$=512$点）
  \item フレーム周期：\qty{10}{m\sec}（$=160$点）
  \item FFT分析点数：512点
  \item 窓の種類：ハニング窓
\end{description}

\section{使用器具}
表\ref{tab:siyoukigu}に使用器具を示す．
\begin{table}[h]
  \centering
  \caption{使用器具}
  \label{tab:siyoukigu}
  \begin{tabular}{c|c|c|c} \hline
    器具名 & 規格 & 製造会社 & その他 \\ \hline 
    パソコン & GALLERIA XF & Dospara &  \\
    OS & Windows 11 25H2 26200.8246 & Microsoft & \\ 
    Cコンパイラ & gcc 8.1.0 & GNU & \\
    グラフ描画ツール & gnuplot 6.0.4 & gnuplot & \\ \hline 
  \end{tabular}
\end{table}

\section{実験結果}
\subsection{バイナリファイルの読み込み}
まず，先頭30ミリ秒が何バイトまでなのかを計算する．
sample01.swは標本化周波数が\qty{16}{k\hertz}であるから，1サンプルあたりの時間は，
\[
  \frac{1}{\qty{16}{k\hertz}} = \qty{0.0625}{\milli\sec}
\]
となる．
先頭30ミリ秒のデータが欲しいから，そのサンプル数は，
\[
  \frac{\qty{30}{\milli\sec}}{\qty{0.0625}{\milli\sec}} = 480
\]
となる．
したがって，480サンプル分のデータを使えば良い．

図\ref{fig:kadai1}にgnuplotを用いて指定のデータを可視化したグラフを示す．
\begin{figure}[h]
  \centering
  \includegraphics[scale=0.4]{kadai1.png}
  \caption{sample01.swの先頭30msecのデータの波形} \label{fig:kadai1}
\end{figure}
6ミリ秒から24ミリ秒にかけて規則的な波形が見られる．
最大振幅の絶対値が1となるように描画されており，正しく出力できたことがわかる．

波形の可視化に必要なデータは，リスト\ref{code:c_kadai}のプログラム（144〜161行目）で生成している．
プログラム内では\texttt{sample01.sw}をバイナリモードで読み込み，先頭から480サンプル分を取得している（153〜159行目）．
また，読み込んだ\qty{16}{bit}整数のデータ値を最大値32768で割ることで，$-1.0$から$1.0$の範囲に正規化した振幅値を求め（157行目），テキストファイル\texttt{kadai1.txt}に出力している．

可視化の際には，リスト\ref{code:gnu_plot}のスクリプトを一括で使用した\cite{gnuplot}．
4〜10行目が課題1のグラフを描画するための記述である．
横軸と縦軸の範囲や名称を設定し，\texttt{plot}コマンドで生成したテキストファイルからデータを読み込み，\texttt{with lines}で点ではなく線でつないで描画している．

\subsection{$DFT$および$FFT$による振幅スペクトルの計算と可視化\cite{fft}}
データ処理を一括で行うプログラム（kadai.c）の実装をリスト\ref{code:c_kadai}に示す．
フーリエ変換における複素数を扱うため，実部\texttt{r}と虚部\texttt{i}をメンバに持つ構造体\texttt{Complex}を定義した（10〜13行目）．
また，C言語標準では複素数演算が直接定義されていないため，以下の基本操作を関数として実装した（15〜32行目）．
\begin{itemize}
    \item \texttt{getComplex}：複素数の生成
    \item \texttt{addComplex}, \texttt{subComplex}, \texttt{timeComplex}：複素数同士の加減乗算
    \item \texttt{divComplexByReal}：複素数の実数による除算
    \item \texttt{cexptheta}：オイラーの公式 $e^{i\theta} = \cos\theta + i\sin\theta$
\end{itemize}

関数\texttt{dft}（88〜97行目）では，離散フーリエ変換の定義式\eqref{eq:dft}を直接計算している．
アルゴリズムの性質上，シグマループと信号長による反復（91〜95行目）という二重ループの計算が必要であり，計算量は$O(N^2)$となる．

関数\texttt{fft}（61〜86行目）では，時間間引き型のCooley-Tukeyアルゴリズムを実装した\cite{kyoto-CT}．
本実装はビット反転とバタフライ演算の2つのステップで構成される．

FFTをインプレース（同一配列内）で計算するため，入力データの添字をビット反転順に並べ替える必要がある．
関数\texttt{bit\_reversal}（42〜59行目）は，信号長$N$が2のべき乗であることを前提とし，ビットシフト演算を用いて高速に並べ替えを実行する．

並べ替えられたデータに対し，バタフライ演算を行う．
外側のループ（63行目〜）でステージ数（$\log_2 N$回）を管理し，内側のループ（69行目〜）で各ブロックの複素積と加減算（75〜76行目）を行う．
回転子$W$の更新には，各ステップで基本となる回転係数\texttt{wm}を順次乗算していく手法（78行目）をとり，三角関数の呼び出し回数を最小限に抑えている．
このアルゴリズムにより，計算量は$O(N \log N)$に削減される．

DFTとFFTの計算効率を定量的に評価するため，Windows環境の提供する高分解能タイマー（\texttt{QueryPerformanceCounter}）を用いて処理時間を計測した\cite{qpc}．
音声データsample01.swの先頭512点をバイナリモードで読み込んで配列の初期値に代入した後，DFT（176〜178行目）およびFFT（181〜183行目）の各アルゴリズムの実行直前と直後のカウント差から秒単位の経過時間を算出した．

このプログラムを使用するには，以下のコマンドを実行する．
\[\texttt{gcc kadai.c -o kadai ; .\textbackslash kadai.exe}\]

計測時間の出力結果をリスト\ref{list:elapsedTime}に示す．
\newpage
\begin{lstlisting}[caption=512点データに対するDFTおよびFFTの実行速度の差（実行出力）,label=list:elapsedTime]
Kadai 2: DFT Elapsed Time: 0.024472 s
Kadai 2: FFT Elapsed Time: 0.000071 s
\end{lstlisting}
FFTに比べDFTは約345倍遅い．

作成したプログラムは，各課題のグラフ描画に必要なすべての振幅スペクトルをテキストファイルとして自動的に出力する仕様とした．
また出力時には，最大振幅値を基準としたdB計算を行う\texttt{outputSpectrum}関数（100〜124行目）を呼び出している．

得られたデータをgnuplotを用いて可視化する．
可視化したグラフを図\ref{fig:kadai2}に示す．
\begin{figure}[tbp]
  \begin{minipage}[b]{0.48\columnwidth}
    \centering
    \includegraphics[width=\columnwidth]{kadai2_linear.png}
    \subcaption{振幅スペクトル（単位なし）}
  \end{minipage}
  \hspace{0.04\columnwidth}
  \begin{minipage}[b]{0.48\columnwidth}
    \centering
    \includegraphics[width=\columnwidth]{kadai2_db.png}
    \subcaption{振幅スペクトル[$\mathrm{dB}$]}
  \end{minipage}
  \caption{sample01.swの振幅スペクトル}\label{fig:kadai2}
\end{figure}
可視化の際には，リスト\ref{code:gnu_plot}のスクリプト（12〜30行目）を使用した．
dB変換の計算は，プログラムの\texttt{outputSpectrum}関数（100〜124行目）で事前に済ませているため，スクリプトではそのまま列を指定してプロットしている．
図を見れば，DFTとFFTのどちらで離散フーリエ変換を行っても，結果は同じになっていることがわかる．

\subsection{窓関数の導入}\label{sec:window}
離散フーリエ変換を行う前の入力信号配列に窓関数をかければ適用できる．
配列に窓関数を適用する関数\texttt{applyWindow}（34〜40行目）を実装した．
配列の全要素に対し，式\eqref{eq:hanning}を計算してかける．
この関数を課題3の処理（193行目）で呼び出し実行している．
得られたデータをgnuplotで可視化した結果を図\ref{fig:kadai2_hanning}に示す．
\begin{figure}[tbp]
  \begin{minipage}[b]{0.48\columnwidth}
    \centering
    \includegraphics[width=\columnwidth]{kadai3_linear.png}
    \subcaption{窓関数適用前後の振幅スペクトル比較（単位なし）}\label{fig:kadai3}
  \end{minipage}
  \hspace{0.04\columnwidth}
  \begin{minipage}[b]{0.48\columnwidth}
    \centering
    \includegraphics[width=\columnwidth]{kadai3_db.png}
    \subcaption{窓関数適用前後の振幅スペクトル比較[$\mathrm{dB}$]}\label{fig:kadai3w}
  \end{minipage}
  \caption{窓関数適用で生じる振幅スペクトルの違い}\label{fig:kadai2_hanning}
\end{figure}
可視化の際には，リスト\ref{code:gnu_plot}のスクリプト（32〜50行目）を使用した．
前述の課題と同様に，dB変換の計算はプログラム内で済ませている．

図\ref{fig:kadai3}および図\ref{fig:kadai3w}を見ると，窓関数を適用したことで，振幅スペクトル全体が小さくなり，特定の周波数のピークがより明瞭になっていることが確認できる．


\subsection{各自の音声の周波数分析}
レポート作成者本人の音声として，「あ」，「い」，「う」，「え」，「お」と発声しているの5つの音声データを用意した．
それぞれ1つの音声データにまとめられているため，各音声の開始時刻を確認し，それに対応する場所から512点を読み込み，窓関数を適用してから離散フーリエ変換を行う．
プログラム（198〜231行目）にて，配列\texttt{targets}に各母音の開始時刻を設定し，ループ処理内で開始時刻からバイト位置を計算している．
217行目で\texttt{fseek}関数を用いて，ファイルポインタを求めたバイト位置まで移動させ，そこを基準に音声データを読み込んでいる．

プログラムを実行し，得られたデータをgnuplotで可視化した結果を図\ref{fig:kadai4}に示す．
\begin{figure}[tbp]
  \begin{minipage}[b]{0.48\columnwidth}
    \centering
    \includegraphics[width=\columnwidth]{kadai4_a_db.png}
    \subcaption{「あ」と発声している音声データ}
  \end{minipage}
  \hspace{0.04\columnwidth}
  \begin{minipage}[b]{0.48\columnwidth}
    \centering
    \includegraphics[width=\columnwidth]{kadai4_i_db.png}
    \subcaption{「い」と発声している音声データ}
  \end{minipage}\\
  \begin{minipage}[b]{0.48\columnwidth}
    \centering
    \includegraphics[width=\columnwidth]{kadai4_u_db.png}
    \subcaption{「う」と発声している音声データ}
  \end{minipage}
  \hspace{0.04\columnwidth}
  \begin{minipage}[b]{0.48\columnwidth}
    \centering
    \includegraphics[width=\columnwidth]{kadai4_e_db.png}
    \subcaption{「え」と発声している音声データ}
  \end{minipage}\\
  \begin{minipage}[b]{\columnwidth}
    \centering
    \includegraphics[width=0.48\columnwidth]{kadai4_o_db.png}
    \subcaption{「お」と発声している音声データ}
  \end{minipage}
  \caption{各音声データの振幅スペクトル[$\mathrm{dB}$]}\label{fig:kadai4}
\end{figure}
可視化の際には，リスト\ref{code:gnu_plot}のスクリプト（52〜72行目）を使用した．
こちらもdB値の計算はプログラム側で行っている．
各振幅スペクトルを見ると，次のことがわかる．
\begin{itemize}
  \item 「あ」：\qty{1}{k\hertz}付近にピークがあり，高域になるにつれて減衰している．
  \item 「い」：\qty{300}{\hertz}付近と，2,3,5\si{k\hertz}付近の高い位置にピークがある．
  \item 「う」：全体的に低い周波数(\qty{2}{k\hertz})以下に集中している．
  \item 「え」：\qty{500}{\hertz}付近にピークがあり，2～3\si{k\hertz}に山がある．
  \item 「お」：\qty{500}{\hertz}と\qty{2.5}{k\hertz}付近にピークがあり，それ以降は落ち込んでいる．
\end{itemize}

\subsection{各自の音声のスペクトログラムの可視化}
音声データ（\texttt{j22330.sw}）の全区間（0〜4.4秒）を対象とし，フレーム長512点，フレーム周期160点ごとにハニング窓を適用してFFTを行い，スペクトログラムを描画した．
図\ref{fig:hatten1}に生成したモノクロおよびカラーのスペクトログラムを示す．
可視化にはリスト\ref{code:gnu_plot}のスクリプト（74〜92行目）を使用し，グローバルな最大振幅を基準（0\si{dB}）として可視化した．

\begin{figure}[h]
  \centering
  \begin{minipage}[b]{0.8\columnwidth}
    \centering
    \includegraphics[width=\columnwidth]{hatten1_mono.png}
    \subcaption{モノクロ表示（黒：強，白：弱）}
  \end{minipage}\\
  \vspace{0.5cm}
  \begin{minipage}[b]{0.8\columnwidth}
    \centering
    \includegraphics[width=\columnwidth]{hatten1_color.png}
    \subcaption{カラー表示（Jet風カラーマップ）}
  \end{minipage}
  \caption{\texttt{j22330.sw}全体のスペクトログラム}\label{fig:hatten1}
\end{figure}

図\ref{fig:hatten1}から，0.3秒〜2.6秒にかけて，"あ","い","う","え","お"の5つの単独発声が，無音区間を挟んで明確に分離して観測される．
それぞれの母音において，低い周波数から高い周波数にかけて特徴的な黒い横帯（フォルマント）が現れており，その高さや間隔が発声内容によって異なっていることが確認できる．
例えば，「あ」は\qty{1}{k\hertz}付近の低い帯と\qty{3}{k\hertz}付近の帯が目立ち，「い」は\qty{300}{\hertz}付近の低い帯と\qty{3}{k\hertz}付近の高い帯の2極に分かれているなど，課題4で得られた振幅スペクトルのピーク位置と一致している．

3.2秒〜4.1秒にかけては，「あいうえお」と連続して発声している区間である．
この連続発声区間では，無音で途切れることなく，時間経過とともにフォルマントの周波数が連続的かつ滑らかに上下に推移していく様子が非常に明確に捉えられている．
これは声道が母音の発声に合わせて滑らかに変化していることを示しており，単独のフレームの周波数分析では得られない，スペクトログラム特有の動的な周波数特性の可視化結果である．

\section{考察}
\subsection{窓関数適用による振幅値の減少}
\ref{sec:window}節で確認した通り，窓関数を適用した音声データの振幅スペクトルは，そのままの信号の振幅スペクトルに比べ全体的に振幅値が小さくなっていることがわかった．
これは窓関数適用によるスペクトル漏れの改善とトレードオフの関係である．
窓関数は式\eqref{eq:hanning}からわかるように，基本的に1以下の値である．
この影響で信号の大部分が減衰させられるため，振幅スペクトルが減少する．

信号のパワーは振幅の2乗で求められる\cite{power}．
音声データの窓関数適用前後で各要素のパワーを足し合わせ，サンプル数で割ることで平均を求め検証する．
検証の計算は，プログラム内の\texttt{calcAmpPower}関数（127〜134行目）を用いて行っている（233〜254行目）．
実行した結果の出力はリスト\ref{code:kosatsu1}のようになった．

\begin{lstlisting}[caption=窓関数適用による信号の平均パワーの変化（実行出力）,label=code:kosatsu1]
Kosatsu 1: BEFORE : 224097506.908203
Kosatsu 1: AFTER  :  60690556.767161
Kosatsu 1: DIFF   : 163406950.141042
Kosatsu 1: RATIO  :         0.270822
\end{lstlisting}
窓関数を適用すると，信号の平均パワーが約0.27倍になることがわかった．

長さ$N$の配列信号$x[n]$の信号の平均パワー$P_x$は以下の式で計算できる．
\[
P_x=\frac{1}{N}\sum_{n=0}^{N-1}\abs{x[n]}^2
\]
次に，この信号に窓関数$w[n]$を適用した信号$y[n]$は，
\[y[n]=x[n]\cdot w[n]\]
となる．
窓関数適用後の信号パワー$P_y$は，
\[
P_y=\frac{1}{N}\sum_{n=0}^{N-1}\abs{x[n]w[n]}^2=\frac{1}{N}\sum_{n=0}^{N-1}\abs{x[n]}^2\abs{w[n]}^2
\]
ここで，ハニング窓関数の定義式\eqref{eq:hanning}より，$w[n]$は
\[
w[n]=0.5-0.5\cos\left(\frac{2\pi n}{N-1}\right)
\]
である．
この関数の値域に関して，$-1 \leqq \cos(\theta) \leqq 1$であり，
\begin{description}
  \item $\cos(\theta)=1$のとき：$w[n]=0.5-0.5\cdot1=0$
  \item $\cos(\theta)=-1$のとき：$w[n]=0.5-0.5\cdot-1=1$
\end{description}
となるから，すべての$n$に対して
\[0\leqq w[n]\leqq1\]
が成り立つ．
これより，$0\leqq \abs{w[n]}^2\leqq1$であるから，以下の関係が成り立つ．
\[\abs{x[n]}^2\abs{w[n]}^2\leqq\abs{x[n]}^2\]
これを全要素について総和を取ると，
\[
\sum_{n=0}^{N-1}\abs{x[n]}^2\abs{w[n]}^2\leqq\sum_{n=0}^{N-1}\abs{x[n]}^2
\]
両辺を$N$で割れば，
\[P_y\leqq P_x\]
となる．

\subsection{ダウンサンプリングの副作用}
より低い標本化周波数に再標本化することをダウンサンプリングという\cite{down}．
\qty{16}{k\hertz}で標本化された信号を，\qty{8}{k\hertz}にダウンサンプリングする方法として「何も前処理をせずに，単純に観測信号を$1/2$に間引く（2回に1回サンプルを破棄する）」がある．
この方法により，高周波成分がスキップされてしまうことによる高音域情報の損失や，信号に標本化周波数の半分より高い周波数成分が含まれているときに起こる折り返し雑音\cite{orikaeshi}といった副作用が生じてしまう．
折り返し雑音が生じると，高周波成分が低周波成分として混入してしまい，本来の音声とは異なるノイズが発生してしまう．
これは標本化定理により説明できる．
例えば，信号に\qty{6}{k\hertz}の正弦波が含まれているとする．
このとき，\qty{16}{k\hertz}で標本化されていれば，$\qty{6}{k\hertz}<\qty{8}{k\hertz}=\frac{\qty{16}{k\hertz}}{2}$であり，標本化定理を満たすから，正しく音を再現することができる．
しかし，ダウンサンプリングを行い，\qty{8}{k\hertz}の信号として再生すると，\qty{2}{k\hertz}の音として観測される．
これが折り返し雑音である．

信号の標本化を数式で表すと，時間$t$は$\frac{n}{f_s}$となるから，元の信号$x[n]$は，
\[x[n]=\cos(2\pi ft)=\cos(2\pi\cdot6000\cdot\frac{n}{16000})=\cos(\frac{3\pi}{4}n)\]
と表される．
$1/2$に間引く操作は$n$を$2m$と置くことで表せるから，ダウンサンプリング後の信号$y[n]$は，
\[y[m]=x[2m]=\cos(\frac{3\pi}{4}\cdot2m)=\cos(\frac{3\pi}{2}m)\]
となる．
ここで，三角関数の性質$\cos(\theta+2\pi n)=\cos(\theta)$と$\cos(-\theta)=\cos(\theta)$を利用すれば，
\[\cos(\frac{3\pi}{2}m)=\cos(2\pi m-\frac{\pi}{2}m)=\cos(-\frac{\pi}{2}m)=\cos(\frac{\pi}{2}m)\]
と変形できる．
この結果を新しい標本化周波数$f_s'=\qty{8000}{\hertz}$で表せば，
\[\cos(2\pi\cdot f_{new}\cdot\frac{n}{8000})=\cos(\frac{\pi}{2}m)\]
より，$f_{new}=\qty{2}{k\hertz}$と求められる．

これらのダウンサンプリングによる副作用を低減させるには，前処理としてローパスフィルタを適用し，標本化周波数の1/2以上の成分を除去する必要がある．
このとき使われるローパスフィルタを，アンチエイリアシングフィルタという\cite{anti}．

\subsection{窓関数と振幅スペクトルの関係\cite{window}}
窓関数の形状は，その後の周波数分析における振幅スペクトルの概形に影響する．
本実験では，窓関数としてハニング窓を扱った．
他にも有名な窓関数を以下に示す．
\begin{itemize}
  \item 矩形窓：窓関数を掛けない状態である．
  \item ハミング窓：両端が完全に0にならない．高い周波数分解能が求められる場合に使う．\[w(x)=0.54-0.46\cos(2\pi x)\]
  \item ブラックマン窓：ハニング窓よりも裾野の広がりが小さい．振幅の小さい信号を大きい信号の近くで測定したい場合に使う．\[w(x)=0.42-0.5\cos(2\pi x)+0.08\cos(4\pi x)\]
\end{itemize}
窓関数は，時間領域で信号との積により適用する．
時間領域の積は周波数領域の畳み込みに対応するから，振幅スペクトルは元の信号のスペクトルと窓関数のスペクトルの畳み込みで得られる．
信号のスペクトルにピークがあれば，窓関数のスペクトルにより畳み込まれ，そのピークが横に広がったようなグラフになる．

窓関数のスペクトルは，スペクトルの中心にある最も大きな山であるメインローブの幅と，その左右に現れる小さな山であるサイドローブのレベルから特徴づけられる．
メインローブの幅は周波数分解能に関わる．
幅が広いほど近接する周波数成分がつながってしまい，区別できなくなる．
サイドローブのレベルはスペクトル漏れに関わる．
サイドローブが高ければ，強い信号の近くにある微小な信号が埋もれて観測できなくなる．

矩形窓は，メインローブが鋭く，サイドローブが高いスペクトルである．
つまり，周波数分解能は高いが，スペクトル漏れが多くなると考えられる．
ハニング窓やハミング窓，ブラックマン窓は，メインローブが太いが，サイドローブが低いスペクトルである．
つまり，周波数分解能が低いが，スペクトル漏れは少ないと考えられる．

\section{まとめ}
本実験では，1時限時系列信号である音声データの分析を対象に，周波数解析の実装，窓関数の影響について理解することができた．
C言語を用いて離散フーリエ変換と高速フーリエ変換を実装し，その計算効率を比較した．
得られたデータをgnuplotを用いて波形観測を行うことができた．
スペクトル漏れを防ぐために窓関数を適用することを知った．
また，その窓関数にも様々な種類が存在して，行うデータ分析によって使用する窓関数を選択する必要があることもわかった．
母音の音声データの周波数分析を行い，発声内容によって振幅スペクトルのピーク位置が異なることを確認した．
これにより，音声の情報の差異が周波数成分の分布として表現されていることを実データから学んだ．


\section{参考文献}
\begin{thebibliography}{9}
\vspace{-1cm}
  \bibitem{gnuplot} Thomas Williams \& Colin Kelley「gnuplot 6.0」，\\ \url{http://www.gnuplot.info/docs/gnuplot-ja.pdf}（更新:2023/12）（参照:2026/04/29）
  \bibitem{fft} tatsukisaito「フーリエ変換を用いて多項式の掛け算を行う」，Teckbot!，\\ \url{https://blog.applibot.co.jp/2020/12/16/fourier-transform/}（更新:2020/12/16）（参照:2026/04/29）
  \bibitem{kyoto-CT} Takuya OOURA「FFT (高速フーリエ・コサイン・サイン変換) の概略と設計法」，京都大学数理解析研究所，\\ \url{https://www.kurims.kyoto-u.ac.jp/~ooura/fftman/}（参照:2026/04/29）
  \bibitem{qpc} 参考書庫「C++ 高速化 処理時間の計測」，\\ \url{http://www.sanko-shoko.net/note.php?id=rnfd}（更新:2018/07/13）（参照:2026/04/29）
  \bibitem{power} MATLABヘルプセンター「信号のパワーの測定」，\\ \url{https://jp.mathworks.com/help/signal/ug/measure-the-power-of-a-signal.html}（参照:2026/04/29）
  \bibitem{down} e-Words「ダウンサンプリング」，\\ \url{https://e-words.jp/w/ダウンサンプリング.html}（参照:2026/04/29）
  \bibitem{orikaeshi} 代表高瀬「ノイズ対策講座-8　折り返し雑音とは？」，デルタテックラボラトリ，\\ \url{https://deltatech-labo.com/2023/04/16/noise_seminar_08}（更新:2025/03/07）（参照:2026/04/29）
  \bibitem{anti} ANALOG DEVICES「フィルタの基礎：アンチエイリアシング」，\\ \url{https://www.analog.com/jp/resources/technical-articles/guide-to-antialiasing-filter-basics.html}（更新:2002/01/11）（参照:2026/04/29）
  \bibitem{window} kennzoの備忘録「窓関数」，\\ \url{https://www.kennzo.net/window-function}（更新:2024/12/06）（参照:2026/04/29）
\end{thebibliography}

\newpage
\begin{lstlisting}[caption=作成した一括gnuplotスクリプト(1/3),label=code:gnu_plot]
set terminal pngcairo size 800,600 font "Arial,12"

# 課題1
set output 'kadai1.png'
set xlabel "Time [msec]"
set ylabel "Normalized Amplitude"
set xrange [0:30]
set yrange [-1:1]
set grid
plot 'kadai1.txt' using 1:2 with lines title "sample01.sw"

# 課題2 (リニア)
set output 'kadai2_linear.png'
set xlabel "Frequency [Hz]"
set ylabel "Amplitude"
set xrange [0:8000]
set autoscale y
set grid
plot 'kadai2_dft.txt' using 1:2 with lines title "DFT", \
     'kadai2_fft.txt' using 1:2 with lines title "FFT"

# 課題2 (dB)
set output 'kadai2_db.png'
set xlabel "Frequency [Hz]"
set ylabel "Amplitude [dB]"
set xrange [0:8000]
set yrange [-90:0]
set grid
plot 'kadai2_dft.txt' using 1:3 with lines title "DFT", \
     'kadai2_fft.txt' using 1:3 with lines title "FFT"

# 課題3 (リニア)
set output 'kadai3_linear.png'
set xlabel "Frequency [Hz]"
set ylabel "Amplitude"
set xrange [0:8000]
set autoscale y
set grid
plot 'kadai2_fft.txt' using 1:2 with lines title "No Window", \
     'kadai3_window.txt' using 1:2 with lines title "Hanning Window"

# 課題3 (dB)
set output 'kadai3_db.png'
set xlabel "Frequency [Hz]"
\end{lstlisting}
\newpage
\addtocounter{lstlisting}{-1}
\begin{lstlisting}[caption=作成した一括gnuplotスクリプト(2/3), firstnumber=45]
set ylabel "Amplitude [dB]"
set xrange [0:8000]
set yrange [-90:0]
set grid
plot 'kadai2_fft.txt' using 1:3 with lines title "No Window", \
     'kadai3_window.txt' using 1:3 with lines title "Hanning Window"

# 課題4 (dB - 各母音個別)
set xlabel "Frequency [Hz]"
set ylabel "Amplitude [dB]"
set xrange [0:8000]
set yrange [-90:0]
set grid

set output 'kadai4_a_db.png'
plot 'kadai4_a.txt' using 1:3 with lines title "/a/"

set output 'kadai4_i_db.png'
plot 'kadai4_i.txt' using 1:3 with lines title "/i/"

set output 'kadai4_u_db.png'
plot 'kadai4_u.txt' using 1:3 with lines title "/u/"

set output 'kadai4_e_db.png'
plot 'kadai4_e.txt' using 1:3 with lines title "/e/"

set output 'kadai4_o_db.png'
plot 'kadai4_o.txt' using 1:3 with lines title "/o/"

# 発展課題1 (スペクトログラム)
set pm3d map
set xlabel "Time [sec]"
set ylabel "Frequency [Hz]"
set xrange [0:*]
set yrange [0:8000]
set zrange [-90:0]
set cbrange [-90:0]
set cblabel "Amplitude [dB]"

# モノクロ表示（大きいところが黒、小さいところが白）
set palette defined (-90 "white", 0 "black")
set output 'hatten1_mono.png'
\end{lstlisting}
\newpage
\addtocounter{lstlisting}{-1}
\begin{lstlisting}[caption=作成した一括gnuplotスクリプト(3/3), firstnumber=87]
splot 'spectrogram.txt' using 1:2:3 notitle

# カラー表示（Jet風カラーマップ）
set palette rgbformulae 22,13,-31
set output 'hatten1_color.png'
splot 'spectrogram.txt' using 1:2:3 notitle
\end{lstlisting}

\begin{lstlisting}[caption=作成したデータ処理プログラム(1/8),label=code:c_kadai]
#define _USE_MATH_DEFINES

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <windows.h>

#define LENGTH 512 /* 信号長 ただし，2の整数乗 */

typedef struct {
    double r; /* 実部 */
    double i; /* 虚部 */
} Complex;  /* 複素数 */

Complex getComplex(double r, double i) {
    return (Complex){r, i};
}
Complex addComplex(Complex a, Complex b) {
    return getComplex(a.r + b.r, a.i + b.i);
}
Complex subComplex(Complex a, Complex b) {
    return getComplex(a.r - b.r, a.i - b.i);
}
Complex timeComplex(Complex a, Complex b) {
    return getComplex(a.r * b.r - a.i * b.i, a.r * b.i + a.i * b.r);
}
Complex divComplexByReal(Complex z, double div) {
    return getComplex(z.r / div, z.i / div);
}
Complex cexptheta(double theta) {
    return getComplex(cos(theta), sin(theta));
}

void applyWindow(Complex *x, int N) {
    for (int i = 0; i < N; i++) {
\end{lstlisting}
\newpage
\addtocounter{lstlisting}{-1}
\begin{lstlisting}[caption=作成したデータ処理プログラム(2/8), firstnumber=36]
        double w = 0.5 - 0.5 * cos(2.0 * M_PI * i / (N - 1));
        x[i].r *= w;
        x[i].i *= w;
    }
}

void bit_reversal(Complex *x, int N) {
    int i, j, k;
    Complex temp;
    j = 0;
    for (i = 0; i < N; i++) {
        if (i < j) {
            temp = x[i];
            x[i] = x[j];
            x[j] = temp;
        }
        k = N >> 1;
        while (k <= j && k > 0) {
            j -= k;
            k >>= 1;
        }
        j += k;
    }
}

void fft(Complex *x, int N, int inverse) {
    bit_reversal(x, N);
    for (int stage = 1; (1 << stage) <= N; stage++) {
        int block = 1 << stage;
        int block2 = block >> 1;
        double theta = (inverse ? 2.0 : -2.0) * M_PI / block;
        Complex wm = cexptheta(theta);

        for (int k = 0; k < N; k += block) {
            Complex w = getComplex(1.0, 0);
            for (int j = 0; j < block2; j++) {
                Complex u = x[k + j];
                Complex v = timeComplex(w, x[k + j + block2]);

                x[k + j] = addComplex(u, v);
                x[k + j + block2] = subComplex(u, v);

                w = timeComplex(w, wm);
            }
\end{lstlisting}
\newpage
\addtocounter{lstlisting}{-1}
\begin{lstlisting}[caption=作成したデータ処理プログラム(3/8), firstnumber=80]
        }
    }
    if (inverse) {
        for (int i = 0; i < N; i++)
            x[i] = divComplexByReal(x[i], N);
    }
}

void dft(Complex *x, Complex *c, int N) {
    for (int k = 0; k < N; k++) {
        c[k] = getComplex(0.0, 0.0);
        for (int n = 0; n < N; n++) {
            double theta = -2.0 * M_PI * k * n / N;
            Complex w = cexptheta(theta);
            c[k] = addComplex(c[k], timeComplex(x[n], w));
        }
    }
}

/* 振幅スペクトルの最大値を0dBとしてテキストファイルに出力する関数 */
void outputSpectrum(const char *fname, Complex *x, int N, double fs) {
    double max_amp = 0.0;
    double *amps = malloc((N / 2 + 1) * sizeof(double));

    // 全データの中で最大振幅を探す
    for (int i = 0; i <= N / 2; i++) {
        amps[i] = sqrt(x[i].r * x[i].r + x[i].i * x[i].i);
        if (amps[i] > max_amp)
            max_amp = amps[i];
    }

    FILE *fp = fopen(fname, "w");
    if (!fp) {
        free(amps);
        return;
    }
    for (int i = 0; i <= N / 2; i++) {
        double f = (double)i * fs / N;
        // 最大振幅を基準(X_b)としてdB値を計算
        double db = (amps[i] > 0 && max_amp > 0) ? 20.0 * log10(amps[i] / max_amp) : -100.0;
        fprintf(fp, "%f %f %f\n", f, amps[i], db);
    }
    fclose(fp);
\end{lstlisting}
\newpage
\addtocounter{lstlisting}{-1}
\begin{lstlisting}[caption=作成したデータ処理プログラム(4/8), firstnumber=123]
    free(amps);
}

/* 複素数配列の平均パワーを求める関数 */
double calcAmpPower(Complex *x, int N) {
    double amp, power = 0.0;
    for (int i = 0; i < N; i++) {
        amp = x[i].r * x[i].r + x[i].i * x[i].i;
        power += amp;
    }
    return power / N;
}

int main(void) {
    int N = LENGTH;
    short input[LENGTH];
    Complex x_dft[LENGTH], x_fft[LENGTH], c[LENGTH];
    FILE *fp;

    double fs = 16000.0;

    printf("Processing Kadai 1...\n");
    fp = fopen("sample01.sw", "rb");
    if (!fp) {
        fprintf(stderr, "Cannot open sample01.sw\n");
        return 1;
    }

    FILE *out1 = fopen("kadai1.txt", "w");
    short val;
    for (int i = 0; i < 480; i++) {
        if (fread(&val, sizeof(short), 1, fp) != 1)
           break;
        double time_ms = (double)i / fs * 1000.0;
        double norm_val = (double)val / 32768.0;
        fprintf(out1, "%f %f\n", time_ms, norm_val);
    }
    fclose(out1);
    fclose(fp);

    printf("Processing Kadai 2...\n");
\end{lstlisting}
\newpage
\addtocounter{lstlisting}{-1}
\begin{lstlisting}[caption=作成したデータ処理プログラム(5/8), firstnumber=164]
    fp = fopen("sample01.sw", "rb");
    fread(input, sizeof(short), N, fp);
    fclose(fp);

    for (int i = 0; i < N; i++) {
        x_dft[i] = getComplex(input[i], 0);
        x_fft[i] = getComplex(input[i], 0);
    }

    LARGE_INTEGER start, end, freq;
    QueryPerformanceFrequency(&freq);

    QueryPerformanceCounter(&start);
    dft(x_dft, c, N);
    QueryPerformanceCounter(&end);
    printf("Kadai 2: DFT Elapsed Time: %f s\n", (double)(end.QuadPart - start.QuadPart) / freq.QuadPart);

    QueryPerformanceCounter(&start);
    fft(x_fft, N, 0);
    QueryPerformanceCounter(&end);
    printf("Kadai 2: FFT Elapsed Time: %f s\n", (double)(end.QuadPart - start.QuadPart) / freq.QuadPart);

    outputSpectrum("kadai2_dft.txt", c, N, fs);
    outputSpectrum("kadai2_fft.txt", x_fft, N, fs);

    printf("Processing Kadai 3...\n");
    Complex x_win[LENGTH];
    for (int i = 0; i < N; i++)
        x_win[i] = getComplex(input[i], 0);
    applyWindow(x_win, N);
    fft(x_win, N, 0);

    outputSpectrum("kadai3_window.txt", x_win, N, fs);

    printf("Processing Kadai 4...\n");
    typedef struct {
        double time_sec;
        const char *vowel;
    } SeekTarget;
\end{lstlisting}
\newpage
\addtocounter{lstlisting}{-1}
\begin{lstlisting}[caption=作成したデータ処理プログラム(6/8), firstnumber=203]
    // 各母音の音声の開始時刻（秒）
    SeekTarget targets[] = {
        {0.361, "a"}, {0.969, "i"}, {1.415, "u"}, {1.936, "e"}, {2.431, "o"}
    };

    fp = fopen("j22330.sw", "rb");
    if (!fp) {
        fprintf(stderr, "Cannot open j22330.sw\n");
        return 1;
    }

    for (int t = 0; t < 5; t++) {
        // 開始時刻からファイル上のバイト位置を計算 (16kHz, 16bit=2byte)
        long start_byte = (long)(targets[t].time_sec * fs * 2);
        fseek(fp, start_byte, SEEK_SET);
        if (fread(input, sizeof(short), N, fp) != (size_t)N)
            continue;

        Complex x_k4[LENGTH];
        for (int i = 0; i < N; i++)
            x_k4[i] = getComplex(input[i], 0);
        applyWindow(x_k4, N);
        fft(x_k4, N, 0);

        char fname[256];
        sprintf(fname, "kadai4_%s.txt", targets[t].vowel);
        outputSpectrum(fname, x_k4, N, fs);
    }
    fclose(fp);

    printf("Processing Kosatsu 1...\n");
    fp = fopen("sample01.sw", "rb");
    if (!fp) {
        fprintf(stderr, "Cannot open sample01.sw\n");
        return 1;
    }
    fread(input, sizeof(short), N, fp);
    fclose(fp);

    Complex x_power[LENGTH];
    for (int i = 0; i < N; i++)
        x_power[i] = getComplex(input[i], 0);

    double beforePower, afterPower;
\end{lstlisting}
\newpage
\addtocounter{lstlisting}{-1}
\begin{lstlisting}[caption=作成したデータ処理プログラム(7/8), firstnumber=247]
    beforePower = calcAmpPower(x_power, N);
    applyWindow(x_power, N);
    afterPower = calcAmpPower(x_power, N);

    printf("Kosatsu 1: BEFORE : %16f\n", beforePower);
    printf("Kosatsu 1: AFTER  : %16f\n", afterPower);
    printf("Kosatsu 1: DIFF   : %16f\n", beforePower - afterPower);
    printf("Kosatsu 1: RATIO  : %16f\n", afterPower / beforePower);

    printf("Processing Hatten Kadai 1...\n");
    fp = fopen("j22330.sw", "rb");
    if (!fp) {
        fprintf(stderr, "Cannot open j22330.sw\n");
        return 1;
    }
    fseek(fp, 0, SEEK_END);
    long file_size = ftell(fp);
    rewind(fp);
    long total_samples = file_size / sizeof(short);
    short *all_input = malloc(file_size);
    if (!all_input) {
        fprintf(stderr, "Memory allocation failed\n");
        fclose(fp);
        return 1;
    }
    fread(all_input, sizeof(short), total_samples, fp);
    fclose(fp);

    int shift = 160; // フレームをずらす間隔 (10ms)
    double global_max = 0.0;
    Complex x_hatten[LENGTH];

    // 0dBの基準となる全体の最大振幅を探索
    for (long i = 0; i <= total_samples - N; i += shift) {
        for (int j = 0; j < N; j++) {
            x_hatten[j] = getComplex(all_input[i + j], 0);
        }
        applyWindow(x_hatten, N);
        fft(x_hatten, N, 0);

        for (int j = 0; j <= N / 2; j++) {
            double amp = sqrt(x_hatten[j].r * x_hatten[j].r + x_hatten[j].i * x_hatten[j].i);
\end{lstlisting}
\newpage
\addtocounter{lstlisting}{-1}
\begin{lstlisting}[caption=作成したデータ処理プログラム(8/8), firstnumber=289]
            if (amp > global_max) {
                global_max = amp;
            }
        }
    }

    // 各フレームのFFTを実行し、最大振幅を基準としたdB値を出力
    FILE *out_sp = fopen("spectrogram.txt", "w");
    if (out_sp && global_max > 0) {
        for (long i = 0; i <= total_samples - N; i += shift) {
            double time_sec = (i + N / 2.0) / fs;
            for (int j = 0; j < N; j++) {
                x_hatten[j] = getComplex(all_input[i + j], 0);
            }
            applyWindow(x_hatten, N);
            fft(x_hatten, N, 0);

            for (int j = 0; j <= N / 2; j++) {
                double f = (double)j * fs / N;
                // 振幅の計算とdB変換 (-90dBで下限クリッピング)
                double amp = sqrt(x_hatten[j].r * x_hatten[j].r + x_hatten[j].i * x_hatten[j].i);
                double db = (amp > 0) ? 20.0 * log10(amp / global_max) : -100.0;
                if (db < -90.0) db = -90.0;
                fprintf(out_sp, "%f %f %f\n", time_sec, f, db);
            }
            fprintf(out_sp, "\n"); // gnuplotのpm3d（面グラフ）描画用に空行を挿入
        }
        fclose(out_sp);
    }
    free(all_input);

    printf("All Kadai processed.\n");
    return 0;
}
\end{lstlisting}


\end{document}