Создание спектрограммы из .wav с использованием БПФ в java

После исследования и множества проб и ошибок я пришел к выводу, что могу построить спектрограмму, которая, как мне кажется, содержит элементы правильного и неправильного.

1. Сначала я считываю файл .wav в массив байтов и извлекаю только часть данных.

2. Я конвертирую байтовый массив в двойной массив, который принимает среднее значение правого и левого каналов. Я также заметил, что 1 выборка 1 канала состоит из 2 байтов. Итак, 4 байта в 1 двойное.

3. Для определенного размера окна степени 2 я применяю БПФ из здесь и получите амплитуду в частотной области. Это вертикальная полоса изображения спектрограммы.

4. Я делаю это несколько раз с тем же размером окна и перекрытием для всех данных и получаю спектрограмму.

Ниже приведен код для чтения .wav в двойной массив.

import java.io.IOException;
import java.nio.ByteBuffer;
import java.nio.file.Files;
import java.nio.file.Path;
import java.nio.file.Paths;
import java.util.Arrays;

public class readWAV2Array {

    private byte[] entireFileData;

    //SR = sampling rate
    public double getSR(){
        ByteBuffer wrapped = ByteBuffer.wrap(Arrays.copyOfRange(entireFileData, 24, 28)); // big-endian by default
        double SR = wrapped.order(java.nio.ByteOrder.LITTLE_ENDIAN).getInt();
        return SR;
    }

    public readWAV2Array(String filepath, boolean print_info) throws IOException{
        Path path = Paths.get(filepath);
        this.entireFileData = Files.readAllBytes(path);

        if (print_info){

        //extract format
        String format = new String(Arrays.copyOfRange(entireFileData, 8, 12), "UTF-8");

        //extract number of channels
        int noOfChannels = entireFileData[22];
        String noOfChannels_str;
        if (noOfChannels == 2)
            noOfChannels_str = "2 (stereo)";
        else if (noOfChannels == 1)
            noOfChannels_str = "1 (mono)";
        else
            noOfChannels_str = noOfChannels + "(more than 2 channels)";

        //extract sampling rate (SR)
        int SR = (int) this.getSR();

        //extract Bit Per Second (BPS/Bit depth)
        int BPS = entireFileData[34];

        System.out.println("---------------------------------------------------");
        System.out.println("File path:          " + filepath);
        System.out.println("File format:        " + format);
        System.out.println("Number of channels: " + noOfChannels_str);
        System.out.println("Sampling rate:      " + SR);
        System.out.println("Bit depth:          " + BPS);
        System.out.println("---------------------------------------------------");

        }
    }

    public double[] getByteArray (){
        byte[] data_raw = Arrays.copyOfRange(entireFileData, 44, entireFileData.length);
        int totalLength = data_raw.length;

        //declare double array for mono
        int new_length = totalLength/4;
        double[] data_mono = new double[new_length];

        double left, right;
        for (int i = 0; i < new_length; i++){
            left = ((data_raw[i] & 0xff) << 8) | (data_raw[i+1] & 0xff);
            right = ((data_raw[i+2] & 0xff) << 8) | (data_raw[i+3] & 0xff);
            data_mono[i] = (left+right)/2.0;
        }       
        return data_mono;
    }
}

Следующий код - это основная программа для запуска

import java.awt.Color;
import java.awt.image.BufferedImage;
import java.io.File;
import java.io.IOException;
import java.util.Arrays;

import javax.imageio.ImageIO;

public class App {

    public static Color getColor(double power) {
        double H = power * 0.4; // Hue (note 0.4 = Green, see huge chart below)
        double S = 1.0; // Saturation
        double B = 1.0; // Brightness

        return Color.getHSBColor((float)H, (float)S, (float)B);
    }

    public static void main(String[] args) {
        // TODO Auto-generated method stub
        String filepath = "audio_work/Sine_Sweep_Full_Spectrum_20_Hz_20_kHz_audiocheck.wav";
        try {

            //get raw double array containing .WAV data
            readWAV2Array audioTest = new readWAV2Array(filepath, true);
            double[] rawData = audioTest.getByteArray();
            int length = rawData.length;

            //initialize parameters for FFT
            int WS = 2048; //WS = window size
            int OF = 8;    //OF = overlap factor
            int windowStep = WS/OF;

            //calculate FFT parameters
            double SR = audioTest.getSR();
            double time_resolution = WS/SR;
            double frequency_resolution = SR/WS;
            double highest_detectable_frequency = SR/2.0;
            double lowest_detectable_frequency = 5.0*SR/WS;

            System.out.println("time_resolution:              " + time_resolution*1000 + " ms");
            System.out.println("frequency_resolution:         " + frequency_resolution + " Hz");
            System.out.println("highest_detectable_frequency: " + highest_detectable_frequency + " Hz");
            System.out.println("lowest_detectable_frequency:  " + lowest_detectable_frequency + " Hz");

            //initialize plotData array
            int nX = (length-WS)/windowStep;
            int nY = WS;
            double[][] plotData = new double[nX][nY]; 

            //apply FFT and find MAX and MIN amplitudes

            double maxAmp = Double.MIN_VALUE;
            double minAmp = Double.MAX_VALUE;

            double amp_square;

            double[] inputImag = new double[length];

            for (int i = 0; i < nX; i++){
                Arrays.fill(inputImag, 0.0);
                double[] WS_array = FFT.fft(Arrays.copyOfRange(rawData, i*windowStep, i*windowStep+WS), inputImag, true);
                for (int j = 0; j < nY; j++){
                    amp_square = (WS_array[2*j]*WS_array[2*j]) + (WS_array[2*j+1]*WS_array[2*j+1]);
                    if (amp_square == 0.0){
                        plotData[i][j] = amp_square;
                    }
                    else{
                        plotData[i][j] = 10 * Math.log10(amp_square);
                    }

                    //find MAX and MIN amplitude
                    if (plotData[i][j] > maxAmp)
                        maxAmp = plotData[i][j];
                    else if (plotData[i][j] < minAmp)
                        minAmp = plotData[i][j];

                }
            }

            System.out.println("---------------------------------------------------");
            System.out.println("Maximum amplitude: " + maxAmp);
            System.out.println("Minimum amplitude: " + minAmp);
            System.out.println("---------------------------------------------------");

            //Normalization
            double diff = maxAmp - minAmp;
            for (int i = 0; i < nX; i++){
                for (int j = 0; j < nY; j++){
                    plotData[i][j] = (plotData[i][j]-minAmp)/diff;
                }
            }

            //plot image
            BufferedImage theImage = new BufferedImage(nX, nY, BufferedImage.TYPE_INT_RGB);
            double ratio;
            for(int x = 0; x<nX; x++){
                for(int y = 0; y<nY; y++){
                    ratio = plotData[x][y];

                    //theImage.setRGB(x, y, new Color(red, green, 0).getRGB());
                    Color newColor = getColor(1.0-ratio);
                    theImage.setRGB(x, y, newColor.getRGB());
                }
            }
            File outputfile = new File("saved.png");
            ImageIO.write(theImage, "png", outputfile);

        } catch (IOException e) {
            // TODO Auto-generated catch block
            e.printStackTrace();
        }
    }

}

Однако изображение, которое я получаю из .wav, воспроизводящего звук с широким диапазоном частот 20–20 кГц, выглядит следующим образом:

Цвет показывает интенсивность звука: красный (высокий) -> зеленый (низкий)

введите описание изображения здесь

По праву, это должно выглядеть примерно так:

введите описание изображения здесь

Я был бы очень признателен, если бы смог получить какое-либо исправление / улучшение / предложение по моему проекту. Заранее благодарю за комментарий по моему вопросу.


person Aung    schedule 02.09.2016    source источник
comment
Можете ли вы объяснить изображение спектрограммы вашей программы: какая ось - время, а какая - частота? Цвета сложно интерпретировать, можно ли использовать оттенки серого?   -  person Ahmed Fasih    schedule 02.09.2016
comment
Ось X - время, а ось Y - частота. красный - высокая амплитуда, зеленый - низкая амплитуда.   -  person Aung    schedule 03.09.2016


Ответы (1)


К счастью, у вас больше прав, чем ошибок.

Первая и основная проблема, которая приводит к появлению лишних красных линий, связана с тем, как вы декодируете данные в readWAV2Array.getByteArray. Поскольку образцы охватывают 4 байта, вы должны индексировать их кратно 4 (например, байты 0,1,2,3 для образца 0, байты 4,5,6,7 для образца 1), иначе вы будете читать перекрывающиеся блоки по 4 байта. (например, байты 0,1,2,3 для образца 0, байты 1,2,3,4 для образца 1). Другая вещь, связанная с этим преобразованием, заключается в том, что вы должны явно привести результат к подписанному типу short, прежде чем его можно будет присвоить left и right (которые относятся к типу double), чтобы получить подписанный 16-битный результат из беззнаковых байтов. Это должно дать вам цикл преобразования, который выглядит так:

for (int i = 0; 4*i+3 < totalLength; i++){
  left = (short)((data_raw[4*i+1] & 0xff) << 8) | (data_raw[4*i] & 0xff);
  right = (short)((data_raw[4*i+3] & 0xff) << 8) | (data_raw[4*i+2] & 0xff);
  data_mono[i] = (left+right)/2.0;
}       

На этом этапе вы должны начать получать график с четкими линиями, представляющими ваш щебетание 20 Гц-20 кГц:

фиксированное декодирование

Но вы должны заметить, что на самом деле вы получаете 2 строки. Это связано с тем, что для сигнала с действительным знаком частотный спектр имеет эрмитову симметрию. Величина спектра выше частоты Найквиста (половина частоты дискретизации, в данном случае 44100 Гц / 2), таким образом, является избыточным отражением спектра ниже частоты Найквиста. Построение только неизбыточной части ниже частоты Найквиста может быть достигнуто путем изменения определения nY в main на:

int nY = WS/2 + 1;

и даст вам:

только неизбыточный спектр

Почти то, что мы ищем, но развертка с увеличением частоты генерирует фигуру с убывающей линией. Это потому, что ваша индексация делает частоту 0 Гц с индексом 0, который является верхней частью рисунка, и частотой 22050 Гц с индексом nY-1, который является нижней частью рисунка. Чтобы перевернуть цифру и получить более обычные 0 Гц внизу и 22050 Гц вверху, вы можете изменить используемую индексацию:

plotData[i][nY-j-1] = 10 * Math.log10(amp_square);

Теперь у вас должен быть график, который выглядит так, как вы ожидали (хотя и с другой цветовой картой):

0 Гц внизу

Последнее замечание: хотя я понимаю ваше намерение избежать логарифма 0 при преобразовании в децибелы, установка выходного сигнала с линейной амплитудой в этом конкретном случае может привести к неожиданным результатам. Вместо этого я бы выбрал амплитуду порога отсечки для защиты:

// select threshold based on the expected spectrum amplitudes
// e.g. 80dB below your signal's spectrum peak amplitude
double threshold = 1.0;
// limit values and convert to dB
plotData[i][nY-j-1] = 10 * Math.log10(Math.max(amp_square,threshold));
person SleuthEye    schedule 03.09.2016
comment
Спасибо, что потратили время на просмотр моего кода и предоставили такое исчерпывающее и подробное объяснение. Я применил коррекцию, но думаю, что форма, которую я получил, не такая, как у вас. Лишние линии убраны, но у меня получилась прямая линия вверх и вниз. Я думаю, что может быть еще одно исправление, которое вы не добавили в комментарий. Если да, то посоветуйте, пожалуйста. Еще раз большое вам спасибо. - person Aung; 04.09.2016
comment
Это изображение, которое я получил после применения всех изменений. Щелкните здесь, чтобы просмотреть изображение. - person Aung; 04.09.2016
comment
Думаю, у меня есть на это ответ. Мой .wav 20-20k Hz - это моно, а не стерео. Я постараюсь исправить это и выложу результат здесь. Спасибо. - person Aung; 04.09.2016
comment
Для тех, кто читает этот пост позже, если файл моно, просто измените 4 на 2, потому что каждые 2 байта - это 1 образец для моно. Большое спасибо @SleuthEye - person Aung; 04.09.2016