You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
322 lines
9.3 KiB
322 lines
9.3 KiB
#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++) {
|
|
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);
|
|
}
|
|
}
|
|
}
|
|
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);
|
|
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");
|
|
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;
|
|
// 各母音の音声の開始時刻（秒）
|
|
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;
|
|
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);
|
|
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;
|
|
}
|
|
|