[DSP编程] C语言谱减法降噪,IFFT后波形出现尖峰脉冲

[复制链接]
387|5
 楼主 | 2018-1-23 14:16 | 显示全部楼层 |阅读模式
本帖最后由 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);
}

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?注册

x
| 2018-1-26 21:03 | 显示全部楼层
把每一帧数据重叠N/2,进行FFT,看看这个现象消失不
 楼主 | 2018-1-27 15:18 | 显示全部楼层
wolf66362629 发表于 2018-1-26 21:03
把每一帧数据重叠N/2,进行FFT,看看这个现象消失不

现在就是重叠N/2呢,N=256,帧移128,问题就是每隔128个点就出现一个尖峰脉冲
| 2018-1-28 18:44 | 显示全部楼层
取数的时候,从N/4到3/4N?
| 2018-4-24 23:21 | 显示全部楼层
在最后的IFFT后做了去窗处理吗?时域重叠部分做了平均处理吗?
 楼主 | 2018-4-29 21:44 | 显示全部楼层
nethopper 发表于 2018-4-24 23:21
在最后的IFFT后做了去窗处理吗?时域重叠部分做了平均处理吗?

你好,IFFT后我有偿试再除以加窗时的函数,但是效果比这还差,时域重叠部分还要做平均处理的吗?
扫描二维码,随时随地手机跟帖
您需要登录后才可以回帖 登录 | 注册

本版积分规则

快速回复

您需要登录后才可以回帖
登录 | 注册
高级模式
我要创建版块 申请成为版主

论坛热帖

关闭

热门推荐上一条 /4 下一条

分享 快速回复 返回顶部 返回列表