本帖最后由 lyh123456 于 2018-1-23 14:16 编辑
相关代码是在网上找的，音频文件格式是8KHz16位，平台是llinux，帧长256，帧移128，加的是汉明窗，过程是先读音频数据，然后进行分帧加窗，FFT,然后求噪声平均值；接着读含噪声语音文件，同样是分帧加窗，FFT,然后求相位并用幅值与噪声幅值相减，最后IFFT，结果波形中每相隔帧移长度的距离就会出现一个尖峰脉冲，对语音波形造成严重失真。
谱减法缺点中的音乐噪声应该是随机性的，现象应该不是这样的；频谱泄漏的话我加了汉明窗并且将FFT点数扩大到了1024字节（此时尖峰脉冲间隔是512），结果还是一样，在此特向论坛中各高手求教。附件是工程所有的文件，具体代码如下： #include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <unistd.h>
#include <sys/types.h>
#include <fcntl.h>
#include <sys/stat.h>
#define WL 256
#define P 10 // 预测系数
#undef pi
#define pi 3.1415926535897932384626434#define winsize 256
#define tempsize winsize/2
#define buffsize winsize*10
typedef struct{
double real;
double img;
}complex;
complex noise[winsize];
double temp[tempsize];
complex x[winsize];
complex NS[16][256];
complex VS[60][256];
complex W[winsize];
complex W1[winsize];
double noise_abs[winsize];
FILE *vfd;
FILE *newfd;
FILE *nfd;
unsigned long wb=0;
unsigned char exitflag=1;
double alpha=10;
double beta=0.02;
double angle(complex a);
double abs1(complex a);
void hamming(complex hw[]);
void del_hamming(complex hw[]);
void hanning(complex hw[]);
void creatreaultfile(char *newpath,char *voicepath);
void noisedata(char *path);
void divi(complex a,complex b,complex *c);
void conjugate_complex(int n,complex in[],complex out[]);
void c_plus(complex a,complex b,complex *c);//复数加
void c_mul(complex a,complex b,complex *c) ;//复数乘
void c_sub(complex a,complex b,complex *c);//复数减法
void c_div(complex a,complex b,complex *c);//复数除法
void fft2(int N,complex f[],char flag);//傅立叶变换 输出也存在数组f中
void ifft2(int N,complex f[]); // 傅里叶逆变换
void c_abs(complex f[],float out[],int n);//复数数组取模
void conjugate_complex(int n,complex in[],complex out[])
{
//complex temp;
int i = 0;
for(i=0;i<n;i++)
{
out.img = -in.img;
out.real = in.real;
}
}
void c_abs(complex f[],float out[],int n)
{
int i = 0;
float t;
for(i=0;i<n;i++)
{
t = f.real * f.real + f.img * f.img;
out = sqrt(t);
}
}
void c_plus(complex a,complex b,complex *c)
{
c->real = a.real + b.real;
c->img = a.img + b.img;
}
void c_sub(complex a,complex b,complex *c)
{
c->real = a.real - b.real;
c->img = a.img - b.img;
}
void c_mul(complex a,complex b,complex *c)
{
c->real = a.real * b.real - a.img * b.img;
c->img = a.real * b.img + a.img * b.real;
}
void c_div(complex a,complex b,complex *c)
{
c->real = (a.real * b.real + a.img * b.img)/(b.real * b.real +b.img * b.img);
c->img = (a.img * b.real - a.real * b.img)/(b.real * b.real +b.img * b.img);
}
#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
void Wn_i(int n,int i,complex *Wn,char flag)
{
Wn->real = cos(2*pi*i/n);
if(flag == 1)
Wn->img = -sin(2*pi*i/n);
else if(flag == 0)
Wn->img = sin(2*pi*i/n);
}
//傅里叶变化
void fft(int N,complex f[],char flag)
{
complex t,wn;//中间变量
int i,j,k,m,n,l,r,M;
int la,lb,lc;
int tt=0;
/*----计算分解的级数M=log2(N)----*/
for(i=N,M=1;(i=i/2)!=1;M++);
//printf("M:%d\n",M);
/*----按照倒位序重新排列原信号----*/
for(i=1,j=N/2;i<=N-2;i++)
{
if(i<j)
{
t=f[j];
f[j]=f;
f=t;
}
k=N/2;
while(k<=j)
{
j=j-k;
k=k/2;
}
j=j+k;
}
/*----FFT算法----*/
for(m=1;m<=M;m++)
{
la=pow(2,m); //la=2^m代表第m级每个分组所含节点数
lb=la/2; //lb代表第m级每个分组所含碟形单元数
//同时它也表示每个碟形单元上下节点之间的距离
/*----碟形运算----*/
for(l=1;l<=lb;l++)
{
r=(l-1)*pow(2,M-m);
for(n=l-1;n<N-1;n=n+la) //遍历每个分组，分组总数为N/la
{
lc=n+lb; //n,lc分别代表一个碟形单元的上、下节点编号
Wn_i(N,r,&wn,flag);//wn=Wnr
c_mul(f[lc],wn,&t);//t = f[lc] * wn复数运算
c_sub(f[n],t,&(f[lc]));//f[lc] = f[n] - f[lc] * Wnr
c_plus(f[n],t,&(f[n]));//f[n] = f[n] + f[lc] * Wnr
tt++;
}
}
}
//printf("tt:%d\n",tt);
}
//傅里叶逆变换
void ifft(int N,complex f[])
{
int i=0;
conjugate_complex(N,f,f);
fft(N,f,1);
conjugate_complex(N,f,f);
for(i=0;i<N;i++)
{
f.img = (f.img)/N;
f.real = (f.real)/N;
}
}
void hamming(complex hw[])
{
double x;
int i;
for(i=0;i<WL;i++)
{
//double cos(x);
x=2*pi*i/(WL-1);
hw.real=(hw.real)*(0.54-0.46*cos(x));//*32768;
}
}
void del_hamming(complex hw[])
{
double x;
int i;
for(i=0;i<WL;i++)
{
x=2*pi*i/(WL-1);
hw.real=(hw.real)/(0.54-0.46*cos(x));
}
}
void hanning(complex hw[])
{
double x;
int i;
for(i=0;i<WL;i++)
{
//double cos(x);
x=2*pi*i/(WL-1);
hw.real=(hw.real)*(0.5-0.5*cos(x));//*32768;
}
}
void divi(complex a,complex b,complex *c)
{
c->real=(a.real*b.real+a.img*b.img)/(b.real*b.real+b.img*b.img);
c->img=(a.img*b.real-a.real*b.img)/(b.real*b.real+b.img*b.img);
}
/*变址计算，将x(n)码位倒置*/
void change(complex *x,int size_x)
{
complex temp;
unsigned short i=0,j=0,k=0;
double t;
for(i=0;i<size_x;i++)
{
k=i;j=0;
t=(log(size_x)/log(2));
while( (t--)>0 )
{
j=j<<1;
j|=(k & 1);
k=k>>1;
}
if(j>i)
{
temp=x;
x=x[j];
x[j]=temp;
}
}
}
double angle(complex a)
{
double m;
//m=atan(a.img)/(a.real);
m=atan(a.img/(a.real+0.000001));
return m;
}
double abs1(complex a)
{
double t;
t = (a.real)*(a.real)+(a.img)*(a.img);
return (double)sqrt(t);
//return t;
}
void readnoisedata(char *path)
{
int i,j,k;
int n1,n2,step=128;
int frame=256;
short buff[2048];
double temp[2048];
nfd=fopen(path,"rb");
if(nfd==NULL)
{
printf("open noise file error!\n");
}
//printf("nfd:%d\n",nfd);
if(fseek(nfd,44,SEEK_SET)==-1)
{
perror("lseek noise file error!\n");
}
if(fread(buff,sizeof(short),2048,nfd)==EOF)
{
perror("read noise file error!\n");
}
for(i=0;i<15;i++)//分帧
{
int n1,n2;
n1=i*step;//i*512
n2=i*step+frame;//i*512+1024
for(j=n1;j<n2;j++)
{
NS[j-n1].real=buff[j];
NS[j-n1].img=0.0;
}
}
fclose(nfd);
}
void creatreaultfile(char *newpath,char *voicepath)
{
char buf[44];
vfd=fopen(voicepath,"r+b");
if(vfd==NULL)printf("open voice file error\n");
fseek(vfd,0,SEEK_SET);
fread(buf,sizeof(char),44,vfd);
newfd=fopen(newpath,"wb");
if(newfd==NULL)
{
printf("creat new file error\n");
}
else
{
if(fwrite(buf,sizeof(char),44,newfd)==EOF)
{
printf("write new file head error\n");
}
}
}
void updateresultfilehead(unsigned long len)
{
unsigned long temptotallen,reallen;
unsigned char headbuf[4];
temptotallen=len-8+44;
headbuf[0]=temptotallen&0xff;
headbuf[1]=temptotallen>>8;
headbuf[2]=temptotallen>>16;
headbuf[3]=temptotallen>>24;
fseek(newfd,4,SEEK_SET);
fwrite(headbuf,sizeof(char),4,newfd);
reallen=len;
headbuf[0]=reallen&0xff;
headbuf[1]=reallen>>8;
headbuf[2]=reallen>>16;
headbuf[3]=reallen>>24;
fseek(newfd,40,SEEK_SET);
fwrite(headbuf,sizeof(char),4,newfd);
fclose(newfd);
}
/**************************主程序*****************************/
void main(void)
{
short tempbuf[7680];
double wavin[7680],wavout[7680];
complex voice[256];
complex voice1[256];
complex noise1[256];
double noise_foward15frame[16][256];
double am_noise[256];
int fs=8000,printfflag=1;
double voice_timedomain[60][256];
double phase[256];
double am_signal[256];
double am_voice[256];
int frame_len=256,step_len=128,n_frame=15,wav_length=7680,i,j,size_x=256,nifrm,ifrm;
int kk=0,n=1;
creatreaultfile("/home/am335x/test/nr/result.wav","/home/am335x/test/nr/test8k.wav");
printf("creat file OK\n");
readnoisedata("/home/am335x/test/nr/noise8k.wav");//read noise data
printf("read noise file OK\n");
switch (fs)
{
case 8000:
frame_len=256;step_len=128;break;
case 10000:
frame_len=400;step_len=200;break;
case 12000:
frame_len=480;step_len=240;break;
case 16000:
frame_len=640;step_len=320;break;
case 44100:
frame_len=1800;step_len=900;break;
default:
frame_len=1800;step_len=900;break;
}
n_frame=(wav_length-frame_len)/step_len+1;//(7680-256)/128+1=58+1=59
//deal the noise data
for(nifrm=0;nifrm<15;nifrm++)//get the noise
{
hamming(NS[nifrm]);
for(i=0;i<frame_len;i++)//256
{
noise1.real=NS[nifrm].real;
noise1.img=0.0;
}
fft(256,noise1,1);
for( i=0;i<frame_len;i++)//256
{
noise_foward15frame[nifrm]=abs1(noise1); //noise_foward15frame保存前15帧的噪音短时傅立叶变换幅度结果 }
}
for(j=0;j<frame_len;j++)
{
double w=0.0;
for(i=0;i<15;i++)
{
w=w+noise_foward15frame[j];
}
am_noise[j]=w/15;
//printf("%f\n",am_noise[j]);
}
exitflag=1;
while(exitflag)
{
i=fread(tempbuf,sizeof(short),wav_length,vfd);
for(i=0;i<wav_length;i++)
{
wavin=tempbuf;
//if(i==256)printf("in:%f \n",wavin);
}
for(i=0;i<n_frame;i++) //分帧59
{
int n1,n2;
n1=i*step_len;//i*256
n2=i*step_len+frame_len;//i*256+512
for(j=n1;j<n2;j++)
{
VS[j-n1].real=wavin[j];
VS[j-n1].img=0.0;
}
}
for(ifrm=0;ifrm<n_frame;ifrm++)//59
{
hamming(VS[ifrm]);
for( i=0;i<frame_len;i++)
{
voice1.real=VS[ifrm].real;
voice1.img=0.0;
}
fft(256,voice1,1);
for(i=0;i<frame_len;i++)
{
phase=angle(voice1); //保存这帧语音信号的傅立叶变换的结果的相
}
for(i=0;i<frame_len;i++)
{
am_signal=abs1(voice1);//保存这帧语音信号的傅立叶变换的结果的幅度
}
for( i=0;i<frame_len;i++)
{
//am_voice=am_signal;
am_voice=am_signal-alpha*am_noise;
//am_voice=sqrt(am_signal*am_signal-am_noise*am_noise); //谱减 %用信号的幅度减去噪声的幅度得到纯净语音的幅度
}
for(i=0;i<frame_len;i++)
{
if(!(am_voice>beta*am_noise))
{
am_voice=beta*am_noise;
}
//if(am_voice<0)
//{
// am_voice=0.0;
//}
}
for(i=0;i<frame_len;i++)
{
voice.real=am_voice; //组合相位与幅度得到去噪后的纯净语音信号
voice.img=phase*(voice.real);
}
ifft(256,voice);
for(i=0;i<frame_len;i++)
{
voice_timedomain[ifrm]=voice.real;// 求这帧纯净语音信号的傅立叶反变换的实部
}
//printf("voice_td:%f\n",voice_timedomain[ifrm][10]);
}
//%求出纯净语音信号的真实幅度
for(i=0;i<wav_length;i++)
{
wavout=0.0;
}
for(i=0;i<n_frame;i++)
{
int m1,m2,m3;
m3=0;
m1=i*step_len;
m2=i*step_len+frame_len;
for(j=m1;j<m2;j++)
{
m3=j-m1;
wavout[j]=wavout[j]+voice_timedomain[m3];
}
}
for(i=0;i<wav_length;i++)
{
tempbuf=(short)wavout;//
}
if(fwrite(tempbuf,sizeof(short),wav_length,newfd)==EOF)
{
printf("write result file data error\n");
}
wb+=wav_length*2;
if(wb>=wav_length*38)exitflag=0;
}
updateresultfilehead(wb);
}
|
|
楼主
标题置顶
标题高亮
;){kind=link}