[低周波騒音検知システム]
 
[構成]
Raspberry Pi2
USB電源KSY UB310-0520
ケース
USB無線LANWLI-UC-G301N
USBハブU2H-TZ420SBK
USB HDDHDJ-UT2.0B
USBマイクHS-MC02UBK

 
[機能概要]
  • 我が家で、夜間「グォングォン」という頭とか耳に圧迫感のある低周波騒音があります。消費者庁でエコキュートやエネファームの低周波騒音を調査しているとのことで、一緒かどうか分かりませんが、うちでも客観的なデータをPiで取得してみます。
  • USBマイクで定期的に一定時間の録音をします。低価格のマイクはノイズ低減のためか低周波帯域がカットされていることが多く、周波数帯域30Hz~のものを見つけて使いましたが、もう生産終了しているようです。
  • 録音したwaveファイルをローパスフィルタで低周波音のみのwaveデータへ変換します。
  • この低周波音のみのwaveデータから平均音量を算出することで、低周波音の大きさを記録します。マイクのA/D変換の仕様が分からないので音量は相対値です。

 
[ローパスフィルタ]
  • 消費者庁のエコキュートの調査データは周波数帯域ごとに測定していますが、観測されたのは40Hz付近のようです。サイン波で音を生成して、うちで聞こえる音と比較してみると、大体80Hzのサイン波と同じに聞こえるので、ちょっと違います。
  • ローパスフィルタは以下のような特性としました。以下は、サイン波で生成した各周波数のwaveファイルを、50Hzをエッジ周波数として、シグマを20Hzとしたローパスフィルタに通した結果です。測定ではエッジ周波数90Hzのフィルタも使いました。
    (参考書:音のプログラミング
    LPF
  • ソースコード(サイン波生成プログラム:makewave.cpp)
    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    #include <string.h>
    
    typedef struct
    {
      int fs; /* sampling frequency */
      int bits; /* data bits */
      int length; /* data length */
      double *s; /* data */
    } MONO_PCM;
    
    void mono_wave_write(MONO_PCM *pcm, char *file_name);
    
    int main(int argc, char **argv)
    {
      MONO_PCM pcm1;
      int n, i, loopnum;
      double f;
      char outfile[20]; 
    
      if(argc>1 && strlen(argv[1])<6){
        f = atoi(argv[1]);
        if(f<10 || f>30000){
          f = 100;
        }
      }else{
        f = 100;
      }
      loopnum = 5;
      pcm1.fs = 8000;
      pcm1.bits = 16;
      pcm1.length = 8000*loopnum;
      pcm1.s = new double[pcm1.length];
    
      for (n=0; n < pcm1.length; n++){
        pcm1.s[n] = sin((double)2.0 * M_PI * f * n / pcm1.fs);
      }
    
      sprintf(outfile, "%d.wav", (int)f); 
      
      mono_wave_write(&pcm1, outfile);
      
      delete[] pcm1.s;
      
      return 0;
    }
    
    void mono_wave_write(MONO_PCM *pcm, char *file_name)
    {
      FILE *fp;
      int n;
      char riff_chunk_ID[4];
      long riff_chunk_size;
      char riff_form_type[4];
      char fmt_chunk_ID[4];
      long fmt_chunk_size;
      short fmt_wave_format_type;
      short fmt_channel;
      long fmt_samples_per_sec;
      long fmt_bytes_per_sec;
      short fmt_block_size;
      short fmt_bits_per_sample;
      char data_chunk_ID[4];
      long data_chunk_size;
      short data;
      double s;
      
      riff_chunk_ID[0] = 'R';
      riff_chunk_ID[1] = 'I';
      riff_chunk_ID[2] = 'F';
      riff_chunk_ID[3] = 'F';
      riff_chunk_size = 36 + pcm->length * 2;
      riff_form_type[0] = 'W';
      riff_form_type[1] = 'A';
      riff_form_type[2] = 'V';
      riff_form_type[3] = 'E';
      
      fmt_chunk_ID[0] = 'f';
      fmt_chunk_ID[1] = 'm';
      fmt_chunk_ID[2] = 't';
      fmt_chunk_ID[3] = ' ';
      fmt_chunk_size = 16;
      fmt_wave_format_type = 1;
      fmt_channel = 1;
      fmt_samples_per_sec = pcm->fs;
      fmt_bytes_per_sec = pcm->fs * pcm->bits / 8;
      fmt_block_size = pcm->bits / 8;
      fmt_bits_per_sample = pcm->bits;
      
      data_chunk_ID[0] = 'd';
      data_chunk_ID[1] = 'a';
      data_chunk_ID[2] = 't';
      data_chunk_ID[3] = 'a';
      data_chunk_size = pcm->length * 2;
      
      fp = fopen(file_name, "wb");
      
      fwrite(riff_chunk_ID, 1, 4, fp);
      fwrite(&riff_chunk_size, 4, 1, fp);
      fwrite(riff_form_type, 1, 4, fp);
      fwrite(fmt_chunk_ID, 1, 4, fp);
      fwrite(&fmt_chunk_size, 4, 1, fp);
      fwrite(&fmt_wave_format_type, 2, 1, fp);
      fwrite(&fmt_channel, 2, 1, fp);
      fwrite(&fmt_samples_per_sec, 4, 1, fp);
      fwrite(&fmt_bytes_per_sec, 4, 1, fp);
      fwrite(&fmt_block_size, 2, 1, fp);
      fwrite(&fmt_bits_per_sample, 2, 1, fp);
      fwrite(data_chunk_ID, 1, 4, fp);
      fwrite(&data_chunk_size, 4, 1, fp);
      
      for (n = 0; n < pcm->length; n++)
      {
        s = (pcm->s[n] + 1.0) / 2.0 * 65536.0;
        
        if (s > 65535.0)
        {
          s = 65535.0;
        }
        else if (s < 0.0)
        {
          s = 0.0;
        }
        
        data = (short)(s + 0.5) - 32768;
        fwrite(&data, 2, 1, fp);
      }
      
      fclose(fp);
    }
    
  • ソースコード(ローパスフィルタ:wavefilter.cpp)
    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    #include <string.h>
    
    void Hanning_window(double w[], int N);
    void FIR_LPF(double fe, int J, double b[], double w[]);
    double sinc(double x);
    
    typedef struct
    {
      int fs; /* sampling frequency */
      int bits; /* data bits */
      int length; /* data length */
      double *s; /* data */
    } MONO_PCM;
    
    int mono_wave_read(MONO_PCM *pcm, char *file_name);
    int mono_wave_write(MONO_PCM *pcm, char *file_name);
    
    int main(int argc, char **argv)
    {
      MONO_PCM pcm0, pcm1;
      int n, m, J;
      double fe, delta, *b, *w, average0, average1;
    
      if(argc>1){
        if(mono_wave_read(&pcm0, argv[1])!=0){
          return 0;
        }
      }else{
        return 0;
      }
      
      pcm1.fs = pcm0.fs;
      pcm1.bits = pcm0.bits;
      pcm1.length = pcm0.length;
      pcm1.s = new double[pcm1.length];
      memset(pcm1.s, 0, pcm1.length);
      
      fe = 50.0 / pcm0.fs;
      delta = 20.0 / pcm0.fs;
      
      J = (int)(3.1 / delta + 0.5) - 1;
      if (J % 2 == 1)
      {
        J++;
      }
      
      b = new double[J + 1];
      w = new double[J + 1];
      
      Hanning_window(w, (J + 1));
      
      FIR_LPF(fe, J, b, w);
     
      average0 = 0; 
      average1 = 0; 
      for (n = 0; n < pcm1.length; n++)
      {
        for (m = 0; m <= J; m++)
        {
          if (n - m >= 0)
          {
            pcm1.s[n] += b[m] * pcm0.s[n - m];
          }
        }
        if(pcm1.s[n]<0){
          average1 -= pcm1.s[n];
        }else{
          average1 += pcm1.s[n];
        }
        if(pcm0.s[n]<0){
          average0 -= pcm0.s[n];
        }else{
          average0 += pcm0.s[n];
        }
      }
      average0 = average0 / pcm0.length;
      average1 = average1 / pcm1.length;
    //  printf("Original average=%f\n",average0);
    //  printf("Filter average=%f\n",average1);
      printf("%f\n",average1);
    //  printf("Filter rate=%.3f\%\n",100*average1/average0);
      
    //  mono_wave_write(&pcm1, (char*)"filter.wav");
      
      delete[] pcm0.s;
      delete[] pcm1.s;
      delete[] b;
      delete[] w;
      
      return 0;
    }
    
    int mono_wave_read(MONO_PCM *pcm, char *file_name)
    {
      FILE *fp;
      int n;
      char riff_chunk_ID[4];
      long riff_chunk_size;
      char riff_form_type[4];
      char fmt_chunk_ID[4];
      long fmt_chunk_size;
      short fmt_wave_format_type;
      short fmt_channel;
      long fmt_samples_per_sec;
      long fmt_bytes_per_sec;
      short fmt_block_size;
      short fmt_bits_per_sample;
      char data_chunk_ID[4];
      long data_chunk_size;
      short data;
      
      fp = fopen(file_name, "rb");
      if(!fp){
        return -1;
      }
      
      fread(riff_chunk_ID, 1, 4, fp);
      fread(&riff_chunk_size, 4, 1, fp);
      fread(riff_form_type, 1, 4, fp);
      fread(fmt_chunk_ID, 1, 4, fp);
      fread(&fmt_chunk_size, 4, 1, fp);
      fread(&fmt_wave_format_type, 2, 1, fp);
      fread(&fmt_channel, 2, 1, fp);
      fread(&fmt_samples_per_sec, 4, 1, fp);
      fread(&fmt_bytes_per_sec, 4, 1, fp);
      fread(&fmt_block_size, 2, 1, fp);
      fread(&fmt_bits_per_sample, 2, 1, fp);
      fread(data_chunk_ID, 1, 4, fp);
      fread(&data_chunk_size, 4, 1, fp);
      
      pcm->fs = fmt_samples_per_sec;
      pcm->bits = fmt_bits_per_sample;
      pcm->length = data_chunk_size / 2;
      pcm->s = new double[pcm->length];
      memset(pcm->s, 0, pcm->length);
      
      for (n = 0; n < pcm->length; n++)
      {
        fread(&data, 2, 1, fp);
        pcm->s[n] = (double)data / 32768.0;
      }
      
      fclose(fp);
      return 0;
    }
    
    int mono_wave_write(MONO_PCM *pcm, char *file_name)
    {
      FILE *fp;
      int n;
      char riff_chunk_ID[4];
      long riff_chunk_size;
      char riff_form_type[4];
      char fmt_chunk_ID[4];
      long fmt_chunk_size;
      short fmt_wave_format_type;
      short fmt_channel;
      long fmt_samples_per_sec;
      long fmt_bytes_per_sec;
      short fmt_block_size;
      short fmt_bits_per_sample;
      char data_chunk_ID[4];
      long data_chunk_size;
      short data;
      double s;
      
      riff_chunk_ID[0] = 'R';
      riff_chunk_ID[1] = 'I';
      riff_chunk_ID[2] = 'F';
      riff_chunk_ID[3] = 'F';
      riff_chunk_size = 36 + pcm->length * 2;
      riff_form_type[0] = 'W';
      riff_form_type[1] = 'A';
      riff_form_type[2] = 'V';
      riff_form_type[3] = 'E';
      
      fmt_chunk_ID[0] = 'f';
      fmt_chunk_ID[1] = 'm';
      fmt_chunk_ID[2] = 't';
      fmt_chunk_ID[3] = ' ';
      fmt_chunk_size = 16;
      fmt_wave_format_type = 1;
      fmt_channel = 1;
      fmt_samples_per_sec = pcm->fs;
      fmt_bytes_per_sec = pcm->fs * pcm->bits / 8;
      fmt_block_size = pcm->bits / 8;
      fmt_bits_per_sample = pcm->bits;
      
      data_chunk_ID[0] = 'd';
      data_chunk_ID[1] = 'a';
      data_chunk_ID[2] = 't';
      data_chunk_ID[3] = 'a';
      data_chunk_size = pcm->length * 2;
      
      fp = fopen(file_name, "wb");
      if(!fp){
        return -1;
      }
      
      fwrite(riff_chunk_ID, 1, 4, fp);
      fwrite(&riff_chunk_size, 4, 1, fp);
      fwrite(riff_form_type, 1, 4, fp);
      fwrite(fmt_chunk_ID, 1, 4, fp);
      fwrite(&fmt_chunk_size, 4, 1, fp);
      fwrite(&fmt_wave_format_type, 2, 1, fp);
      fwrite(&fmt_channel, 2, 1, fp);
      fwrite(&fmt_samples_per_sec, 4, 1, fp);
      fwrite(&fmt_bytes_per_sec, 4, 1, fp);
      fwrite(&fmt_block_size, 2, 1, fp);
      fwrite(&fmt_bits_per_sample, 2, 1, fp);
      fwrite(data_chunk_ID, 1, 4, fp);
      fwrite(&data_chunk_size, 4, 1, fp);
      
      for (n = 0; n < pcm->length; n++)
      {
        s = (pcm->s[n] + 1.0) / 2.0 * 65536.0;
        
        if (s > 65535.0)
        {
          s = 65535.0;
        }
        else if (s < 0.0)
        {
          s = 0.0;
        }
        
        data = (short)(s + 0.5) - 32768;
        fwrite(&data, 2, 1, fp);
      }
      
      fclose(fp);
      return 0;
    }
    
    void FIR_LPF(double fe, int J, double b[], double w[])
    {
      int m;
      int offset;
      
      offset = J / 2;
      for (m = -J / 2; m <= J / 2; m++)
      {
        b[offset + m] = 2.0 * fe * sinc(2.0 * M_PI * fe * m);
      }
      
      for (m = 0; m < J + 1; m++)
      {
        b[m] *= w[m];
      }
    }
    
    double sinc(double x)
    {
      double y;
      
      if (x == 0.0)
      {
        y = 1.0;
      }
      else
      {
        y = sin(x) / x;
      }
      
      return y;
    }
    
    void Hanning_window(double w[], int N)
    {
      int n;
      
      if (N % 2 == 0)
      {
        for (n = 0; n < N; n++)
        {
          w[n] = 0.5 - 0.5 * cos(2.0 * M_PI * n / N);
        }
      }
      else
      {
        for (n = 0; n < N; n++)
        {
          w[n] = 0.5 - 0.5 * cos(2.0 * M_PI * (n + 0.5) / N);
        }
      }
    }
    
    
  • データ取得はcronとperlスクリプトで行います。
    cron
    # /etc/cron.d/wavefilter
    
    # record low freq. sound every 5 minutes
    0-59/5 *     * * *     root   cd /home/pi/cpp/wavefilter; ./run_filter.sh
    
  • perl(run_filter.sh)
    #! /usr/bin/perl
    
    my ($sec, $min, $hour, $mday, $mon, $year, $wday, $yday, $isdst) = localtime(time);
    $wavefile = sprintf("%02d%02d%02d.wav", $hour, $min, $sec);
    $timetxt = sprintf("%04d/%02d/%02d %02d:%02d:%02d,", $year+1900, $mon+1, $mday, $hour, $min, $sec);
    `arecord -r 8000 -f S16_LE -d 5 -D plughw:2 $wavefile`;
    `echo -n "$timetxt" >> /home/pi/cpp/wavefilter/record.log`;
    `./wavefilter $wavefile >> /home/pi/cpp/wavefilter/record.log`;
    exit;
    
[測定結果]
  • 音量は、waveデータの最大値を1として0.0001を基準としたdB値にしています。この絶対値に意味はありません。
    Result
    結果としては、夜間の低周波騒音を目に見える形には出来ませんでした。昼間にランダムに検出している音をフィルタを通して確認してみると、歩いたりドタバタしているときの床の振動音に含まれる低周波成分で、低周波音の検知は出来ているようです。やはり問題の音は音量としては小さく、少なくとも昼間の生活音に含まれる低周波音よりはるかに小さいレベルのようです。騒音の原因は不明ですが、こうなると気にし過ぎと言われても言い返せず。
 
 
TOP PAGE