스펙트로그램의 개요
목소리나 음악 등의 소리를 수학적으로 기술하는 방법으로는, 크기를 나타내는 시그널을 시간에 대한 함수로 구하는 것이 있는데요. 이를 직접 들여다보는 것 보다 주파수 별로 그 크기가 어떻게 달라지는지를 분석하면, 소리의 패턴을 과학이나 공학적으로 이해하기가 더 쉬워집니다. 스펙트로그램 (spectrogram)은 서로 다른 주파수의 기여분이 시간에 따라 어떻게 달라지는지를 이미지로 나타낸 것입니다.
Spectrogram - Wikipedia
From Wikipedia, the free encyclopedia Jump to navigation Jump to search Visual representation of the spectrum of frequencies of a signal as it varies with time Spectrogram of the spoken words "nineteenth century". Frequencies are shown increasing up the ve
en.wikipedia.org
비단 소리뿐만 아니라 모든 종류의 시그널에 대해 스펙트로그램을 만드는 것이 가능합니다만, 주로 음향학이나 언어학 및 광학에서 사용된다고 하는군요. 예를 들어서 패턴분석을 통한 음성인식을 구현하거나, 레이더 기술에도 응용이 된다고 합니다.
시그널의 주파수 같은 개념이 낮설게 느껴지신다면, 파동의 개념에 대한 다음 포스팅을 읽어보면 도움이 될 거라고 생각합니다. 뿐만 아니라 파동의 주기와 주파수, 파장과 파동수 사이의 상관관계 등에 대해서도 포괄적으로 다루었기 때문에, 흔하게 볼 수 있는 자연현상인 파동과 좀 더 친숙해질 수 있는 계기가 되리라 믿습니다.
물리학 상식 : 파동의 이해
이번 포스팅에서는 파동과 이를 구성하는 주파수, 파장, 진폭 등의 기본 개념들을 자세하게 짚어보겠습니다. 파동은 자연계에서 매우 흔하게 발견되는 패턴인데요. 우리가 일상생활에서 가장
swstar.tistory.com
시간에 대한 함수를 주파수에 대한 함수로 변환하는 작업이기 때문에, 수학적으로는 푸리에 변환에 그 기반을 두고 있는데요. 한가지 눈여겨볼 점이 있다면, 스펙트로그램은 시간과 주파수 모두를 인자 (parameter)로 하는 시그널의 크기라는 것입니다. 개인적으로 처음에 좀 헷갈렸던 부분인데, 왜냐하면 푸리에 변환으로 연결되는 두 함수는 시간과 주파수 중의 하나만을 인자로 받기 때문이죠. 시간에 따라 주파수 스펙트럼이 어떻게 변화하는지를 나타내기 위해서, 시그널에다가 윈도우 함수 (window function)라는 것을 별도로 곱한 뒤에 푸리에 변환을 수행한다고 합니다.
이 윈도우 함수가 시간에 따라 증가하다가 특정 시각에 피크가 되고, 그 이후로는 감소하는 형태를 생각해볼 수 있습니다. 예를 들면 가우스 함수가 이러한 형태를 띄고 있죠. 그렇게 되면 윈도우 함수가 최대가 되는 시점에서 시그널의 주파수 스펙트럼을 볼 수 있게 되고, 이를 시각화 한 것이 바로 스펙트로그램입니다.
수학적으로는 단기 푸리에 변환 (Short-time Fourier Transform, 줄여서 STFT)을 수행한 뒤에, 변환된 복소수 함수의 절대값의 제곱으로 시그널의 세기를 구합니다. 일반적으로 가로축은 시각, 즉 단기 푸리에 변환에서 윈도우 함수가 최대가 되는 시각이 됩니다. 세로축은 주파수에 해당되고, 시그널의 세기를 서로 다른 색깔, 명도, 채도 등으로 나타냅니다.
주어진 시그널을 가지고 스펙트로그램을 구하는 것을 간단히 도식화하면 다음과 같습니다.
여기서는 설명을 쉽게 하기 위해서 시간별로 서로 다른 색깔을 사용하고, 시그널의 세기를 채도로 나타냈는데요. 실제로는 다수의 스펙트로그램에서 색깔이나 밝기로 시그널의 세기를 나타냅니다.
C++ 프로그램
이번에는 시간에 대한 함수로 주어지는 시그널을 가지고, 스펙트로그램의 정보를 구하는 C++ 프로그램을 소개해 볼까 합니다. 그에 앞서서 먼저 푸리에 변환을 수행하기 위한 라이브러리가 필요한데, 여기서는 고속 푸리에 변환을 위한 FFTW를 사용하도록 하겠습니다. FFTW 라이브러리의 설치 및 사용법은 다음 포스팅에 소개되어 있습니다.
FFTW : C/C++ 고속 푸리에 변환 라이브러리
라이브러리 소개 현재 진행중인 연구에서 사용중인 C언어 푸리에 변환 라이브러리인 FFTW에 대해 간략히 포스팅해볼까 합니다. FFTW는 Fastest Fourier Transform in the West 의 약자인데요. 저 같은 토종 아
swstar.tistory.com
단기 푸리에 변환을 수행하여 스펙트로그램을 얻기 위한 C++ 프로그램을 소개해 봅니다. 먼저 변환 대상이 되는 함수의 포인터를 매개변수로 받아서 푸리에 변환을 수행하는 FFTransformer 라는 클래스의 헤더 및 소스 파일입니다.
FFTransformer.h [다운로드]
FFTransformer.cpp [다운로드]
C언어 또는 C++ 프로그램에서 함수를 다른 함수의 매개변수로 넘겨주는 법에 대해서 포스팅을 한 적이 있는데요. 여기서도 FFTransformer 클래스를 예시로 소개한 바 있습니다.
C/C++ 에서 함수를 매개변수로 사용하기
함수 포인터를 이용한 구현 일반적으로 C언어나 C++ 에서 사용하는 함수의 경우, 인자(매개변수) 혹은 파라미터로 변수를 받아갑니다. 이 값들을 가지고 정의된 기능을 수행하게 되죠. 하지만 프
swstar.tistory.com
이제 시그널과 윈도우 함수의 정보를 이용해 단기 푸리에 변환을 수행하는 소스코드가 필요합니다.
test1_spectrogram_fftw.cpp [다운로드]
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>
#include<vector>
#include"FFTransformer.h"
/* define a signal multiplied by
* a gaussian window function */
class InfoSignal {
private :
bool have_signal_;
public :
/* number of timestep
* which is same with total number of samples */
int nstep_;
double t_ini_; // initial time
double t_fin_; // final time
double delta_t_; // size of timestep
// array to store the signal
std::vector<double> signal_;
/* mean and width
* of the gaussian window function */
double win_mean_;
double win_sigma_;
InfoSignal() {
win_mean_ = 0.;
win_sigma_ = 1.;
have_signal_ = false;
}
~InfoSignal() {}
void set_signal(double t_ini_in, double t_fin_in,
std::vector<double> &signal_in) {
nstep_ = static_cast<int>(signal_in.size());
if (nstep_ < 1) {
return;
}
t_ini_ = t_ini_in;
t_fin_ = t_fin_in;
delta_t_ =
fabs(t_fin_ - t_ini_) / static_cast<double>(nstep_);
signal_.clear();
for (int it = 0; it < nstep_; it++) {
signal_.push_back(signal_in.at(it));
}
signal_.push_back(signal_in.at(0));
have_signal_ = true;
}
bool have_signal(){return have_signal_;}
};
// gaussian function
double func_gaussian(double x, double mean_x, double sigma_x);
/* the real and imaginary parts of a complex function,
* to be Fourier-transformed */
double func_fft_re(double t, void *ptr_param);
double func_fft_im(double t, void *ptr_param);
// the main function
int main(int argc, char *argv[]) {
InfoSignal info_signal;
// width of the window function
info_signal.win_sigma_ = 1.;
/* range in time
* here from 0 to 100 seconds */
double t_min = 0.;
double t_max = 100.;
/* total number of samples
* here 1000 Hz times 100 seconds */
int n_sample = 100 * 1000;
// array to store the signal
std::vector<double> vec_signal;
/* specify the signal
* here a sum of sinusoidal functions */
vec_signal.clear();
for (int isam = 0; isam < n_sample; isam++) {
double t_now = t_min +
static_cast<double>(isam) * (t_max - t_min) /
static_cast<double>(n_sample);
double signal_now = -0.2 +
0.2 * cos(0.25 * t_now) + 0.3 * sin(0.5 * t_now) +
0.5 * cos(0.75 * t_now) - cos(t_now) +
cos(2. * t_now);
vec_signal.push_back(signal_now);
}
fprintf(stdout,
" 1 t_min = %e, t_max = %e\n", t_min, t_max);
fprintf(stdout,
" vec_signal size = %d\n",
static_cast<int>(vec_signal.size()));
info_signal.set_signal(t_min, t_max, vec_signal);
fprintf(stdout,
" 2 info_signal.nstep_ = %d\n", info_signal.nstep_);
fprintf(stdout,
" info_signal.t_ini_ = %e\n", info_signal.t_ini_);
fprintf(stdout,
" info_signal.t_fin_ = %e\n", info_signal.t_fin_);
fprintf(stdout,
" info_signal.delta_t_ = %e\n", info_signal.delta_t_);
// Initialize the Fourier transformation class
FFTransformer fft_transform;
bool fft_init =
fft_transform.init(info_signal.t_ini_, info_signal.t_fin_,
info_signal.delta_t_, false);
if (!fft_init) {
fprintf(stderr,
"ERROR : FFTransformer is not properly initialized.\n");
exit(1);
}
fprintf(stdout, " 3 FFTransformer initialized\n");
// function pointer
PtrFuncFFT ptr_func_in;
ptr_func_in.re = &func_fft_re; // real part
ptr_func_in.im = &func_fft_im; // imaginary part
/* parameter
* which defines the signal and window function */
ptr_func_in.param = &info_signal;
char filename_spec[200];
char buf_win_sigma[20];
sprintf(buf_win_sigma, "%f", info_signal.win_sigma_);
strcpy(filename_spec, "spectrogram_win_sigma_");
strcat(filename_spec, buf_win_sigma);
strcat(filename_spec, ".dat");
FILE *fout_spec;
fout_spec = fopen(filename_spec, "w");
fprintf(fout_spec,
"# [time] [frequency] [amplitude square]\n");
fprintf(stdout, " 4 performing STFT\n");
int nstep_one_spectrogram = 20;
for (int it = 0; it < info_signal.nstep_; it++) {
double t_now = info_signal.t_ini_ +
info_signal.delta_t_ * static_cast<double>(it);
info_signal.win_mean_ = t_now;
if (it % nstep_one_spectrogram != 0 ||
t_now > 20.) {
continue;
}
// Perform fast Fourier transformation
bool fft_next = fft_transform.next(ptr_func_in);
if (!fft_next) {
fprintf(stderr,
"ERROR : FFTransformer failed to perform.\n");
exit(1);
}
std::vector<double> vec_omega;
std::vector<double> x_stft_re;
std::vector<double> x_stft_im;
// get the result of short-time Fourier transform
fft_transform.get_array_fin(vec_omega,
x_stft_re, x_stft_im);
int n_frequency = static_cast<int>(vec_omega.size());
for (int ifreq = 0; ifreq < n_frequency; ifreq++) {
double freq_now = vec_omega.at(ifreq) / (2. * M_PI);
if (freq_now < 0. || freq_now > 1.) {
continue;
}
double x_re = x_stft_re.at(ifreq);
double x_im = x_stft_im.at(ifreq);
double sqr_amp =
fabs(x_re * x_re) + fabs(x_im * x_im);
fprintf(fout_spec,
" %e %e %e\n", t_now, freq_now, sqr_amp);
}
fprintf(fout_spec, "\n");
}
fprintf(stdout, " 5 finished\n");
fclose(fout_spec);
return 0;
}
double func_gaussian(double x, double mean_x, double sigma_x) {
double dx_scaled = (x - mean_x) / sigma_x;
double func =
exp(-0.5 * fabs(dx_scaled * dx_scaled)) /
(sigma_x * sqrt(2. * M_PI));
return func;
}
/* We are doing Fourier transform
* of the signal multiplied by a gaussian window function. */
double func_fft_re(double t, void *ptr_param) {
InfoSignal *ptr_sig = (InfoSignal *)ptr_param;
if (!ptr_sig->have_signal()) {
return 0.;
}
/* If the time is too far from the peak
* of the window function, we get zero. */
if (fabs(t - ptr_sig->win_mean_) >
5. * ptr_sig->win_sigma_) {
return 0.;
}
int it =
static_cast<int>(floor((t - ptr_sig->t_ini_) /
ptr_sig->delta_t_));
if (it < 0 || it >= ptr_sig->nstep_) {
return 0.;
}
// linear interpolation
double frac_interpolation[2];
frac_interpolation[1] =
(t - ptr_sig->t_ini_) / ptr_sig->delta_t_ -
static_cast<double>(it);
frac_interpolation[0] = 1. - frac_interpolation[1];
double func_signal = 0.;
for (int jt = 0; jt < 2; jt++) {
func_signal +=
frac_interpolation[jt] * ptr_sig->signal_.at(it + jt);
}
double func_window = func_gaussian(t,
ptr_sig->win_mean_, ptr_sig->win_sigma_);
return func_signal * func_window;
}
// And the imaginary part would be zero.
double func_fft_im(double t, void *ptr_param) {
return 0.;
}
시그널 및 윈도우 함수에 대한 정보를 담기 위한 InfoSignal 클래스를 새롭게 도입했습니다. 이 객체는 시그널을 저장한 벡터 및 윈도우 함수의 평균값과 폭 등의 정보를 담고 있죠. 푸리에 변환의 대상이 되는 함수의 실수 (func_fft_re) 및 허수 (func_fft_im) 부분을 정의하는 데 있어서 InfoSignal 객체의 포인터가 들어갑니다. 뿐만 아니라 푸리에 변환을 수행하는 FFTransformer 객체의 인자로 들어가는 포인터 구조체 PtrFuncFFT 가 있는데요. 이것 역시 InfoSignal 객체의 포인터를 가지게 됩니다.
다루고자 하는 시그널에 맞게 초기 시각, 마지막 시각을 정하고 시그널을 벡터에 저장한 다음 InfoSignal 객체의 set_signal 함수를 호출하여 초기화합니다. 첨부된 코드는 아래의 삼각함수 시그널을 위한 것입니다만, 목적에 따라 변형이 가능한 부분입니다.
포인터 구조체 PtrFuncFFT 에 변환대상이 되는 함수와 시그널에 대한 정보를 저장하고, FFTransformer 객체의 init 함수를 호출하여 초기화하면 푸리에 변환을 위한 준비가 끝납니다. 그리고 for 반복문을 통해 각 시간별로 단기 푸리에 변환을 수행하고, 스펙트로그램을 계산합니다.
다음과 같이 소스 파일들을 컴파일하여 실행파일을 얻을 수 있습니다.
g++ FFTransformer.cpp -c
g++ test1_spectrogram_fftw.cpp -c
g++ test1_spectrogram_fftw.o FFTransformer.o -lfftw3 -lm -o [실행파일 이름]
예시 : 삼각함수 시그널
간단한 예시로서 삼각함수들의 합으로 이루어진 시그널을 가지고 스펙트로그램을 한 번 만들어 봅시다. 이번에 다루어 볼 시그널을 그래프로 그려보면 다음과 같습니다.
윈도우 함수의 폭이 1초인 경우를 상정해서 계산을 해 보면, 다음과 같은 스펙트로그램을 얻을 수 있는데요. 위에서 언급한 대로 가로축과 세로축은 각각 시간 및 주파수입니다. 그리고 시그널의 세기에 따라 서로 다른 색깔을 띄고 있습니다.
여기서 한 가지 짚고 넘어갈 점은 윈도우 함수의 폭에 따라 스펙트로그램의 형태가 달라진다는 것입니다. 이를 좀 더 자세히 살펴보기 위해 윈도우 함수의 폭을 0.2초로 줄여서 계산을 해 보면, 다음과 같은 스펙트로그램이 나오게 됩니다.
스펙트로그램이 세로로 길게 뻗어있는 것을 볼 수 있는데요. 이말인즉슨, 시간에 따라 시그널의 세기가 어떻게 변하는지가 잘 드러납니다. 반면에 어느 주파수에서 시그널이 전달되고 있는지는 알기 어렵죠.
이번에는 반대로 윈도우 함수의 폭이 매우 큰 경우를 상정해서 스펙트로그램을 만들어 봅시다. 폭을 5초로 설정하고 계산을 하면 다음과 같은 모양을 얻을 수 있습니다.
윈도우 함수의 폭이 0.2초인 경우와는 반대로, 스펙트로그램이 가로로 길게 뻗어 있는 모양을 하고 있습니다. 그래서 시간에 따라 시그널의 세기가 어떻게 달라지는지를 파악하기가 어렵습니다. 그 대신에 시그널이 가진 주파수들을 알아내기는 상당히 쉽습니다.
이렇게 시간과 주파수 사이의 불확정성이 존재하기 때문에, 하나의 윈도우 함수만 가지고는 시그널의 실체를 제대로 이해할 수 없습니다. 그래서 실무적으로는 여러 개의 다른 폭을 가진 윈도우 함수들로부터 스펙트로그램을 얻은 다음, 이를 취합하여 시그널을 분석한다고 하는군요.
