Skip to content

Instantly share code, notes, and snippets.

@naoh16

naoh16/emsnr.c

Created February 3, 2012 16:57
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save naoh16/1731120 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save naoh16/1731120 to your computer and use it in GitHub Desktop.
Download ZIP
Calculate SNR (Signal-to-Noise Ratio) based on an assumption of two-mixture Gaussian distribution
Raw
emsnr.c
/**
* emsnr.c
*
* SNRを計算する
*
* Usage: echo datafilename | emsnr residual maxloop
* Output: filename SNR low_power high_power low_variance high_variance
*
* HARA Sunao
* 2004-10-31 尤度がNaNの時は終了するように変更
* 2004-10-12 分散の下限を設定
* 2004-09-24 メモリ開放忘れ，ファイル閉じ忘れ修正
* 2004-05-13
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define TPI (double)6.283185307179586476925286766559
#define VAR_FLOOR 1.0E-4
//#define FRAME_SHIFT 100
//#define WINDOW_SIZE 256
#define FRAME_SHIFT 32
#define WINDOW_SIZE 64
// defined in stdio.h
//#define FILENAME_MAX 1024
double g_hamming_table[WINDOW_SIZE];
typedef struct {
double *data;
long len;
double mu1;
double mu2;
double sigma1;
double sigma2;
double lambda;
} GAUSS_PROB_PARAM;
/**
* short(16bit)のデータをdouble型にして取得
*/
double* i16d_readdata(char *filename, long *len, int bSwap) {
short *tmp_data;
double *data, *p;
long i, _len;
// char buf[256];
FILE *fp;
fp = fopen(filename, "rb");
if(!fp){
fprintf(stderr, "Error at fopen(i16d_readdata).\n");
return NULL;
}
// 長さの取得
{
fseek(fp, 0L, SEEK_END);
_len=ftell(fp) / sizeof(short);
*len = _len;
fseek(fp, 0L, SEEK_SET);
}
// RIFF ヘッダのチェックとすっとばし
if(_len>4) {
unsigned char buf[5];
unsigned long chanklen;
fread(buf, sizeof(char), 4, fp); // "RIFF"
buf[4] = '\0';
if(buf[0]=='R' && buf[1]=='I' && buf[2]=='F' && buf[3]=='F') {
fprintf(stderr, "Caution: RIFF File...");
fread(buf, sizeof(char), 4, fp); // "RIFF" length
fread(buf, sizeof(char), 4, fp); // "WAVE"
fread(buf, sizeof(char), 4, fp); // "fmt "
while(!(buf[0]=='d' && buf[1]=='a' && buf[2]=='t' && buf[3]=='a')) {
// printf("%s", buf);
fread(&chanklen, sizeof(long), 1, fp); // "????" length
fseek(fp, chanklen, SEEK_CUR);
fread(buf, sizeof(char), 4, fp); // "????"
}
// printf("%s", buf);
fread(&chanklen, sizeof(long), 1, fp); // "????" length
// printf("[%d->%d]\n", _len, chanklen/sizeof(short));
_len = chanklen / sizeof(short);
*len = _len;
fprintf(stderr, "header was removed.\n");
} else {
// 元に戻しておく
fseek(fp, 0L, SEEK_SET);
}
}
// メモリ確保
tmp_data = malloc(_len * sizeof(short));
if(tmp_data == NULL){
fprintf(stderr, "Error at malloc_tmpdata(i16d_readdata).\n");
return NULL;
}
data = malloc(_len * sizeof(double));
if(data == NULL){
fprintf(stderr, "Error at malloc_data(i16d_readdata).\n");
free(tmp_data);
return NULL;
}
if(bSwap) {
unsigned char *p, tmp;
// 読み込み
fread(tmp_data, sizeof(short), _len, fp);
p = (unsigned char *)tmp_data;
// swap
for(i=0; i<_len; i++){
tmp = p[2*i];
p[2*i] = p[2*i+1];
p[2*i+1] = tmp;
}
} else{
// 読み込み
fread(tmp_data, sizeof(short), _len, fp);
}
// doubleの配列にコピー(代入)
for(i=0; i<_len; i++){
data[i] = (double)tmp_data[i];
}
free(tmp_data); tmp_data = NULL;
fclose(fp);
return data;
}
void init_hamming() {
int n;
double len_1 = (double)(WINDOW_SIZE - 1);
for(n=0; n<WINDOW_SIZE; n++) {
g_hamming_table[n] = (0.54 - 0.46 * cos(TPI*(double)n / len_1));
}
}
double* dd_hamming(double* src, double *dest, long len) {
long n;
// double len_1 = (double)(len - 1);
for(n=0; n<len; n++) {
// dest[n] = src[n] * (0.54 - 0.46 * cos(TPI*(double)n / len_1));
dest[n] = src[n] * g_hamming_table[n];
// printf("%f, ", dest[n]);
}
return dest;
}
double dd_avepower(double* src, long len) {
double result = 0.0;
long i;
for(i=0; i<len; i++){
result += src[i]*src[i] / (double)(len);
}
return result;
}
double* dd_frame_power(double* src, long src_len, long* dest_len) {
long f = 0;
long frame_len;
double *frame_power;
double *frame_tmp;
frame_len = (long)((double)(src_len - WINDOW_SIZE) / (double)FRAME_SHIFT + 0.5);
frame_power = (double *)malloc(frame_len*sizeof(double));
frame_tmp = (double *)malloc(WINDOW_SIZE*sizeof(double));
for(f=0; f<frame_len; f++){
dd_hamming(&src[f*FRAME_SHIFT], frame_tmp, WINDOW_SIZE);
frame_power[f] = dd_avepower(frame_tmp, WINDOW_SIZE);
frame_power[f] = 10.0 * log10(frame_power[f]);
// printf("%d:%f\n", f, frame_power[f]);
}
*dest_len = frame_len;
free(frame_tmp); frame_tmp = NULL;
return frame_power;
}
double gauss_prob(double x, double mu, double sigma){
double x_mu = x-mu;
return exp(-(x_mu*x_mu)/(2*sigma)) / sqrt(TPI*sigma);
}
double mix2_gauss_prob(double x, GAUSS_PROB_PARAM* param){
double x_mu1, x_mu2;
// x_mu1 = x - param->mu1;
// x_mu2 = x - param->mu2;
return param->lambda * gauss_prob(x, param->mu1, param->sigma1)
+ (1-param->lambda) * gauss_prob(x, param->mu2, param->sigma2);
}
void init_gaussparam(GAUSS_PROB_PARAM *param, double* data, long len){
long i;
double mu1, mu2, x_mu;
double sigma1, sigma2;
// 重みの初期値は等分(=0.5)
param->lambda = 0.5;
// 平均の初期値は学習データの最小値と最大値
mu1 = mu2 = data[0];
for(i=1; i<len; i++){
if(data[i] < mu1) mu1 = data[i];
if(data[i] > mu2) mu2 = data[i];
}
param->mu1 = mu1;
param->mu2 = mu2;
// 分散の初期値は1
param->sigma1 = 1.0;
param->sigma2 = 1.0;
}
double loglikelyhood(GAUSS_PROB_PARAM* param, double* data, int data_len){
int i;
double sum=0.0;
for(i=0; i<data_len; i++){
sum += log(mix2_gauss_prob(data[i], param));
}
return sum;
}
double em_2mixture(GAUSS_PROB_PARAM *param, int loop_count, double apply_residual){
int i, k;
int same_count=0;
double new_loglikelyhood = 0.0, old_loglikelyhood = 0.0;
double *data = param->data;
long len = param->len;
// パラメータ初期化
init_gaussparam(param, data, len);
// printf("#loop lambda mu1 sigma1 mu2 sigma2 loglikelyhood\n");
new_loglikelyhood = loglikelyhood(param, data, len);
// printf("%3d %f %f %f %f %f %f\n",
// -1, param.lambda, param.mu1, param.sigma1, param.mu2, param.sigma2, new_loglikelyhood);
for(k=0; k<loop_count; k++){
double prob, tmp_prob;
double prob_sum1=0.0, prob_sum2=0.0;
double new_mu1=0.0, new_mu2=0.0;
double new_sigma1=0.0, new_sigma2=0.0;
double new_lambda=0.0;
double x_mu;
// 重みの最推定
for(i=0; i<len; i++){
tmp_prob = mix2_gauss_prob(data[i], param);
prob = param->lambda * gauss_prob(data[i], param->mu1, param->sigma1);
prob_sum1 += prob / tmp_prob;
prob = (1 - param->lambda) * gauss_prob(data[i], param->mu2, param->sigma2);
prob_sum2 += prob / tmp_prob;
}
param->lambda = prob_sum1 / (prob_sum1 + prob_sum2);
if(param->lambda < 0.00001) param->lambda = 0.00001;
else if(param->lambda > 0.99999) param->lambda = 0.99999;
// 平均の最推定(estimate)
for(i=0; i<len; i++){
tmp_prob = mix2_gauss_prob(data[i], param);
prob = param->lambda * gauss_prob(data[i], param->mu1, param->sigma1);
new_mu1 += data[i] * prob / tmp_prob;
prob = (1 - param->lambda) * gauss_prob(data[i], param->mu2, param->sigma2);
new_mu2 += data[i] * prob / tmp_prob;
}
new_mu1 /= prob_sum1;
new_mu2 /= prob_sum2;
param->mu1 = new_mu1;
param->mu2 = new_mu2;
// 分散の最推定
for(i=0; i<len; i++){
tmp_prob = mix2_gauss_prob(data[i], param);
prob = param->lambda * gauss_prob(data[i], param->mu1, param->sigma1);
x_mu = data[i] - new_mu1;
new_sigma1 += x_mu * x_mu * prob / tmp_prob;
prob = (1 - param->lambda) * gauss_prob(data[i], param->mu2, param->sigma2);
x_mu = data[i] - new_mu2;
new_sigma2 += x_mu * x_mu * prob / tmp_prob;
}
param->sigma1 = new_sigma1 / prob_sum1;
param->sigma2 = new_sigma2 / prob_sum2;
// 分散フロアリング
if(param->sigma1 < VAR_FLOOR) param->sigma1 = VAR_FLOOR;
if(param->sigma2 < VAR_FLOOR) param->sigma2 = VAR_FLOOR;
new_loglikelyhood = loglikelyhood(param, data, len);
#if _DEBUG
printf("%3d %f %f %f %f %f %f\n",
k, param->lambda, param->mu1, param->sigma1, param->mu2, param->sigma2, new_loglikelyhood);
#endif
// NaN 対策
if(isnan(new_loglikelyhood)) break;
if(new_loglikelyhood - old_loglikelyhood < apply_residual){
same_count++;
}else{
same_count = 0;
}
old_loglikelyhood = new_loglikelyhood;
if(same_count > 10) break;
}
return new_loglikelyhood;
}
int main(int argc, char* argv[])
{
int loop_count;
char filename[FILENAME_MAX], *p;
double *data, *data_raw;
double loglikelyhood = 0.0;
double apply_residual = 0.0;
int bSwap = 0;
long len_raw, len;
GAUSS_PROB_PARAM param;
long filecount = 0;
double snrsum = 0.0;
if(argc <= 2){
fprintf(stderr, " usage: echo datafilename | emsnr residual maxloop\n");
fprintf(stderr, "example: echo test.raw | emsnr 1.0E-8 1000\n");
fprintf(stderr, "\nOutput: filename SNR low_power high_power low_variance high_variance loglikelyhood\n");
return 1;
}
// filename = argv[1];
apply_residual = atof(argv[1]);
loop_count = atoi(argv[2]);
if(argc > 3)
bSwap = (strcmp(argv[3], "-x")==0)?1:0;
// ハミング窓テーブル作成
init_hamming();
// データ読み込み
while(NULL != fgets(&filename[0], FILENAME_MAX, stdin)){
double snr = 0.0;
while(NULL != (p=strrchr(filename, '\r'))) *p = '\0';
while(NULL != (p=strrchr(filename, '\n'))) *p = '\0';
printf("%s ", filename);
data_raw = i16d_readdata(filename, &len_raw, bSwap);
if(data_raw == NULL) {
fprintf(stderr, "Error at readdata[%s].\n", filename);
return 1;
}
// frame_powerをdataに
data = dd_frame_power(data_raw, len_raw, &len);
// fprintf(stderr, "frame len: %d\n", len);
printf("%d ", len);
fflush(stdout);
param.data = data; param.len = len;
loglikelyhood = em_2mixture(&param, loop_count, apply_residual);
snr = fabs(param.mu2 - param.mu1);
printf("%f %f %f %f %f %f\n", snr, param.mu1, param.mu2,
param.sigma1, param.sigma2,
loglikelyhood);
fflush(stdout);
free(data_raw); data_raw = NULL;
free(data); data = NULL;
memset(&filename[0], 0, FILENAME_MAX);
// 最後の表示用の保存
if(!isnan(snr)) {
++filecount;
snrsum += snr;
}
}
fprintf(stderr, "Filecount=%d MeanSNR=%.2f\n", filecount, snrsum/filecount);
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment