Фильтр Винера и FFTW

#fft #qt #cpp

Идея реализовать фильтр Винера для восстановления изображений.
Формула здесь.
На сколько понимаю, то необходимо сделать дискретное преобразование фурье ядра и
самого изображения, которое необходимо восстановать, а затем уже с ними работать.
Ниже код, как я думаю, который должен реализовать DFT для входного изображения.
Однако полученные результаты не совпадают с примерами из сети.
Должно было получиться link text, а получилось link text.
//fftw_complex *inputImageFFT;     
//fftw_complex *outputImageFFT;
//fftw_plan plan;  
int N;
width = this->inImg.width();
height = this->inImg.height();
N=width*height;

inputImageFFT = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*N);
outputImageFFT = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*N);
plan = fftw_plan_dft_2d(width,height,inputImageFFT,outputImageFFT,
FFTW_FORWARD,FFTW_ESTIMATE);

   for (int i = 0, k = 0; i < height; i++)
   {
       for (int j = 0; j < width; j++, k++)
       {
           inputImageFFT[k][0] = (uint)inImg.pixel(j,i);
           inputImageFFT[k][1] =  0.0;
        }

   }

   fftw_execute(plan);

   float max=0;
   float mag=0;
   for (int i = 1, k = 1; i < height; i++)
   {
       for (int j = 1; j < width; j++, k++)
       {
            mag = qSqrt(qPow(outputImageFFT[k][0],2) + pow(outputImageFFT[k][1],2));
            if (max < mag)
            max = mag;
       }

  }

   outImg = new QImage(width,height, QImage::Format_RGB32);

   for (int i = 0, k = 0; i < height; i++)
   {
       for (int j = 0; j < width; j++, k++)
       {
           float mag = sqrt(pow(outputImageFFT[k][0],2) + pow(outputImageFFT[k][1],2));
           mag = 255*(mag/max);
           outImg->setPixel(j,i,qRgb(mag,mag,mag));
       }
   }
//вывод outImg

Подскажите, в том ли направлении реализации фильтра я двигаюсь и в какую сторону
рыть?    

Ответы

Ответ 1

Есть отличный набор статей по теме на хабре:
Восстановление расфокусированных и смазанных изображений и
Восстановление расфокусированных и смазанных изображений. Практика  и есть ещё одна
из серии, в статье найдёте ссылку.
Там всё расписано как и что делать. В качестве примера приводится программа, написанная
автором под названием smart deblur.

БПФ звукового сигнала (C++ и SFML)

#cpp #алгоритм #аудио #sfml #fft

У меня вопрос по БПФ звукового сигнала. Я хочу нарисовать график на основе звукового
сигнала. Но я не силен в этой теме. Хочу совета у опытных товарищей с чего начать.

У меня несколько вопросов:
1) Что такое samples (семплы)? Например при использовании мультимедийной библиотеки
SFML для получения семплов используется следующая конструкция:

sf::Buffer Buffer;
Buffer.loadFromFile("sound.wav");
const sf::Int16 *input = Buffer.getSamples();


Так вот я так понимаю семплы это бинарное представление звукового файла? Я правильно
понимаю, что вышепредставленный код, это тоже самое что и:

    typedef short int16;
    int16 *load()
    {
      FILE  *fp;
      if((fp=fopen("sound.wav", "rb"))==NULL) {
        printf("Ошибка при открытии файла.\n");
      }
      fseek(fp, 0, SEEK_END);
      long N = ftell(fp);
      fseek(fp, 0, SEEK_SET);
      int16 *A = new int16[N];
      for(i=0; i

Ответы
Ответ 1

Если речь идет о несжатом звуковом файле, то семпл - это отсчет, полученный при оцифровке
сигнала, т.е. просто мгновенное значение амплитуды аналогового сигнала. По этим отсчетам
можно построить вашу "синусойду", т.е. представление сигнала во временной области.
БПФ же, это быстрое преобразование Фурье, выполнив его мы получаем отображение сигнала
в частотной области (разложение сигнала по частотам). Как реализована конкретная библиотека
я не знаю, но для алгоритма БПФ необходима кратность 2. Хотя можно и нолики в конце
дописать, ну да ладно, это уже ЦОС.
Рекомендую вам для начала ознакомится с теоретической стороной вопроса, чтобы четко
понимать, что вы делаете. Успехов)

Быстрое преобразование Фурье. Как выделить частоту ноты?

#java #аудио #fft #музыка

Заинтересовался реализацией выделением нот и распознаванием сигналов исходя из нот.
Вот если взять фортепиано, записать звук пары клавиш, после чего применить БПФ, на
выходе мы получим массив комплексных чисел, где амплитудой будет модуль комплексного
числа, а аргумент его фазой. Вот немного не понятно как выделить частоту ноты?
    

Ответы

Ответ 1

Преобразование фурье производится на комплексных числах.
На вход для преобразования следует подавать в реальную часть амплитуду(величину сигнала),
а в мнимую часть нуль.

На выходе мы получим массив комплексных чисел, где амплитудой будет модуль комплексного
числа, а аргумент его фазой. Каждый элемент массива представляет из себя одну гармонику,
начиная с нулевой и завершая n-ой. Или это называется спектрами? В разных источниках
на этот счёт разная информация.

Спектры/гармоники отделены друг от друга дискретным шагом равным частотам дискретизации/количество
отсчётов. Количество отсчётов равно количеству чисел амплитуд на входе, длине массивов
входящих комплексных чисел, количеству длины входного сэмпла. В общем, количеству комплексных
чисел на входе.

Больше на входе отсчётов, больше разрешение по частоте на выходе.

То есть, на выходе получив 1024 значений комплексных получаем 1024 значений амплитуда-фаза.
Из этих значений, располагая массив по порядку и амплитуде можно получить нечто вроде
визуализации амплитуды звука по частоте.

Если дискретизация 44100Гц, а входной массив имеет длину 65536 то шаг между элементами
массива на выходе получается 0.672912598 Гц.

В целях распознавания человеческой речи подобная точность бессмысленна и избыточна.
Путём генетического алгоритма (естественного отбора человеческого вида в условиях среды
планеты Земля) оптимальным максимальным будет шаг между частотами от 1 до 9 Гц, то
есть на вход желательно подавать максимум 8192 значений амплитуды-времени для частоты
дискретизации 44100.

Но прежде чем подавать эти данные на распознавание, с ними нужно ещё немного повозиться,
но я не знаю как. Там что-то про окошки и что-то с функциями надо делать. Ничего не понял.

UPD: Короче, вот код уточнения частоты основанный на смещении фаз в дополнение к
фурье на java.

/*числа пи*/
 public static final double SinglePi = Math.PI;
 public static final double DoublePi = 2*Math.PI;

/**
* На вход подаётся два спектра, смещение между ними и частота.
* @param spectrum0
* @param spectrum1
* @param shiftPerFrame
* @param sampleRate
* @return dictionary
*/
public static HashMap GetJoinedSpectrum(
 List spectrum0,
 List spectrum1,
 double shiftsPerFrame,
 double sampleRate)
{
 int frameSize = spectrum0.size();
 double frameTime = frameSize/sampleRate;
 double shiftTime = frameTime/shiftsPerFrame;
 double binToFrequancy = sampleRate/frameSize;
 HashMap dictionary = new HashMap(){
 for (int bin = 0; bin < frameSize; bin++)
 {
  double omegaExpected = DoublePi*(bin*binToFrequancy); // ω=2πf
  double omegaActual = (spectrum1.get(bin).phase() - spectrum0.get(bin).phase())/shiftTime;
// ω=∂φ/∂t
  double omegaDelta = Align(omegaActual - omegaExpected, DoublePi); // Δω=(∂ω + π)%2π - π
  double binDelta = omegaDelta/(DoublePi*binToFrequancy);
  double frequancyActual = (bin + binDelta)*binToFrequancy;
  double magnitude = spectrum1.get(bin).abs()+ spectrum0.get(bin).abs();
  dictionary.put(frequancyActual, magnitude*(0.5 + Math.abs(binDelta)));
 }
 return dictionary;
}

public static double Align(double angle, double period) {
 int qpd = (int) (angle/period);
 if (qpd >= 0) qpd += qpd & 1;
 else qpd -= qpd & 1;
 angle -= period*qpd;
 return angle;
}


Ссылки:

https://m.habrahabr.ru/post/247385/

http://introcs.cs.princeton.edu/java/32class/Complex.java

Алгоритмы активного шумоподавления

#алгоритм #аудио #fft #обработка_сигналов

Меня интересуют алгоритмы активного шумоподавления. Углубляясь в детали реализации,
выполняется спектральное вычитание "на ходу".

Имеется некоторый профиль шума и небольшой длины запись, из которой этот шум нужно
удалить ("вычесть").

Так у меня вопрос по конкретному алгоритму спектрального вычитания: как оно выполняется
и, собственно, над чем? В общих чертах я себе это представляю как вычитание из спектра
(соответствие частот амплитудам гармоник) сигнала спектра шума и обратное преобразование
Фурье. Вообще, вроде как, применяемое здесь FFT даёт не соответствие частот амплитудам,
а просто массив комплексных чисел. Как тут вычитать что из чего не представляю...

P.S.: Не нужно давать ссылки на курсы по обработке звука и спектральному анализу
сигналов. Просто вопрос - просто ясный ответ.



Есть такие специализированные программы, как, например, Noise Gator. Так вот, если
кому известно, подскажите, какие алгоритмы шумоподавления ими используются (в частности,
используется ли спектральное вычитание; если да, то откуда программа берёт слепок шума).
    

Ответы

Ответ 1

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

Любой анализ сигнала с помощью преобразований Фурье предполагает, что сигнал является
стационарным на исследуемом отрезке. Поэтому, если сигнал нестационарный, он разбивается
на отдельные отрезки, называемые окнами. Выбор размера окна зависит от типа исследуемого
сигнала, обычно около 25 - 50 мс для звукового сигнала (меньшие значения - для человеческой
речи, большие - для музыки, особенно состоящей из струнных смычковых инструментов).
Можно использовать перекрывающиеся окна, для повышения точности анализа.

Однако, просто так применить преобразование Фурье к обрезанным окнам нельзя, при
этом некорректно обрабатываются граничные области отрезков. Для решения этой проблемы
сигнал предварительно умножают на специальную весовую функцию ("оконную"). Примеры
оконных функций см. в статье Оконное преобразование Фурье 

Далее выполняется непосредственно преобразование Фурье. Оно дает в результате спектр
сигнала, т.е. значение комплексной амплитуды для различных диапазонов частот. Из него
и надо вычитать спектр шума. Из модуля комплексной амплитуды сигнала вычитается модуль
комплексной амплитуды шума, умноженный на некий коэффициент; если результат отрицательный,
он заменяется на ноль. Фазовый компонент оставляется нетронутым. К результату вычитания
можно применить обратное преобразование Фурье, и получить "очищенный" сигнал.

Итоговый алгоритм шумоподавления при наличии известного образца шума:


Разделение сигнала на окна
Применение оконного преобразования Фурье к окнам
Вычитание (по модулю) спектра амплитуды шума из спектра амплитуды сигнала:



  A = Max ( A с. - k * А ш. ; 0)


где k - коэффициент, подбираемый опытным путем


Применение обратного преобразования Фурье к результату


Размер окна, перекрытие окон, тип применяемой оконной функции подбираются опытным путем.

Ссылки

Removing noise from audio using Fourier transform in Matlab

How can I select an optimal window for Short Time Fourier Transform?


How to select frequency resolution and window size in FFT?


Ответ 2

Вы уже задавали этот вопрос, а я отвечал. Попытаюсь ответить лучше, чем в прошлый раз.


  ОКНО


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



Кажется, что это синус, но на самом деле это вот что:

Резкий перепад, выделенный красным - это очень сильное искажение. Если вы прослушаете
это в наушниках ушами, то поймете, что это совсем не тот шум, который вы хотите вычесть
из своего сигнала. Чтобы уменьшить искажение, применяется окно (не ликвидировать, ведь
окно само вносит искажения в сигнал):



Результат применения окна:



Это была работа с образцом шума.

Теперь можно выполнить БПФ на этом образце шума и получить его спектр. Дальше, на
каждом входном блоке сигнала можно так же выполнить БПФ и получить его спектр.

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

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

Получение частоты звука с микрофона

#c_sharp #аудио #фурье #fft #naudio

Есть задача получать частоту звука с микрофона для дальнейших преобразований на C#.
Подобное уже делал на Python с numpy, но тут как то не клеится...

private void write(byte[] angles, int byte_len)
{
    Complex[] _fftBuffer = new Complex[byte_len];
    var _m = (int)Math.Log(byte_len, 2.0);

    for (var n=0; n < byte_len; n++)
    {
        var r = angles[n];
        var i = 0;
        _fftBuffer[n].X = (float)(r * FastFourierTransform.HammingWindow(n, byte_len));
        _fftBuffer[n].Y = i;

    }
    FastFourierTransform.FFT(true, _m, _fftBuffer);
    float[] fft_x = new float[_fftBuffer.Length];
    for (var i=0; i<_fftBuffer.Length; i++)
    {
        fft_x[i] = Math.Abs(_fftBuffer[i].X);
    }
    int i_peak = fft_x.ToList().IndexOf(fft_x.Max());
    for (var i = 0; i < _fftBuffer.Length; i++)
    {
        fft_x[i] = (float)Math.Log(Math.Abs(_fftBuffer[i].X));
    }
    ///var i_peak = fft_x.Max();
    var i_interp = parabolic(fft_x, i_peak);

    float freq = byte_len * i_interp / angles.Length;
    Console.WriteLine("Debug stop");

    ///this.port.Write(this.first_command, 0, 8);
    ///this.port.Write(this.get_comand(0.0f, 0.0f), 0, 32);


}
private float parabolic(float[] f, int peak)
{
    var xv = 0.5f * (f[peak-1] - f[peak+1])/(f[peak-1]-2*f[peak] + f[peak+1]) + peak;

    return xv;
}


Это что на С# сделал. Сделано так, ибо так же было на Python.
В конечном итоге работает, но неправильно - частота получается одна и та же (+- пара
герц), но по идее там должны быть абсолютно другие частоты - передаем данные с помощью
звука.

Помогите разобраться, что в коде может быть не так?

UPD

Вот что получилось в итоге. Спасибо товарищу @MSDN.WhiteKnight - натолкнул на правильные
мысли.
Плюс использовался проект 
вот отсюда

выкладываю только код, который несет смысл по вытаскиванию частот звука с микрофона(МОНО).
Можно переделать на стерео - не особо сложно будет

public partial class MainWindow : Window
{
    static double Fs = 48000; // Частота дискретизации !В данной программе ТОЛЬКО
целые числа
    static double T = 1.0 / Fs; // Шаг дискретизации
    static int N; //Длина сигнала (точек)
    static double Fn = Fs / 2;// Частота Найквиста
    WaveIn waveIn;

    public MainWindow()
    {
        InitializeComponent();
    }

    void waveIn_DataAvailable(object sender, WaveInEventArgs e)
    {

        //данные из буфера распределяем в массив чтобы в нем они были в формате ?PCM?
        byte[] buffer = e.Buffer;
        N = buffer.Length;
        int bytesRecorded = e.BytesRecorded;
        Complex[] sig = new Complex[bytesRecorded / 2];
        for (int i = 0, j = 0; i < e.BytesRecorded; i += 2, j++)
        {
            short sample = (short)((buffer[i + 1] << 8) | buffer[i + 0]);
            sig[j] = sample / 32768f;
        }

        Fourier.Forward(sig, FourierOptions.Matlab);
        // обнуляем спектр на небольших частотах (там постоянная составляющая и вообще
много помех)
        for (int i = 0; i < 35 * sig.Length / Fn; i++)
        {
            sig[i] = 0;
        }

        write(sig);

    }
    //Окончание записи
    private void waveIn_RecordingStopped(object sender, EventArgs e)
    {
        waveIn.Dispose();
        waveIn = null;
    }

    private void start_button_Click(object sender, RoutedEventArgs e)
    {
        this.waveIn = new WaveIn();
        this.waveIn.DeviceNumber = 0;
        this.waveIn.DataAvailable += this.waveIn_DataAvailable;
        this.waveIn.RecordingStopped += this.waveIn_RecordingStopped;
        this.waveIn.WaveFormat = new WaveFormat((int)Fs, 1);
        this.waveIn.StartRecording();
        Start_button.IsEnabled = false;
        Stop_button.IsEnabled = true;
        this.log_box("старт записи");
    }


    private void stop_Button_Click(object sender, RoutedEventArgs e)
    {
        this.stop_recording();
    }

    private void stop_recording()
    {
        this.waveIn.StopRecording();
        Start_button.IsEnabled = true;
        Stop_button.IsEnabled = false;
        this.log_box("конец записи");
    }


    private void log_box(string message)
    {
        Log_Box.AppendText("\n" + message);
        Log_Box.ScrollToEnd();
    }

    private void write(Complex[] signal)
    {
        PointPairList list1 = new PointPairList();
        int max_index = 0;
        double freq = 0;
        double K = signal.Length / 2;
        for (int i = 0; i < K; i++)
        {
            list1.Add(i * Fn / K, Complex.Abs(signal[i]) / N * 2);
        }

        foreach (ZedGraph.PointPair i in list1)
        {
            if (i.Y > list1[max_index].Y)
            {
                max_index = list1.IndexOf(i);
            }
        }
        freq = list1[max_index].X;

        string s = freq.ToString();
        log_box(s);


    }

}


что использовалось...
NAudio - для получения потока звука с микрофона

MathNET - Фурье

ZedGraph - как переходник для работы с сигналом после преобразования Фурье - его
можно убрать, но в моем случае был удобен.
    

Ответы

Ответ 1

Ваш код будет работать нормально, только если на вход подать данные определенного
формата: моно, 1 байт на сэмпл, определенная частота дискретизации и т.п. Кроме того,
он не учитывает несколько деталей: из результата БПФ нужно отбросить первое значение
("постоянная составляющая") и вторую половину значений (которая не несет полезной информации);
количество сэмплов должно быть в степени 2. 

Лучше написать код, который может корректно обрабатывать разные форматы, для этого
возьмем за основу класс SampleAggregator из примера на Github:

using System;
using System.Collections.Generic;
using System.Text;
using System.Diagnostics;
using NAudio.Dsp;

namespace WindowsFormsTest1
{    
    public class SampleAggregator
    {
        // volume
        public event EventHandler MaximumCalculated;
        private float maxValue;
        private float minValue;
        public int NotificationCount { get; set; }
        public Complex[] FftBuffer { get { return this.fftBuffer; } }
        int count;

        // FFT
        public event EventHandler FftCalculated;
        public bool PerformFFT { get; set; }
        private Complex[] fftBuffer;
        private FftEventArgs fftArgs;
        private int fftPos;
        private int fftLength;
        private int m;

        public SampleAggregator(int fftLength = 1024)
        {
            if (!IsPowerOfTwo(fftLength))
            {
                throw new ArgumentException("FFT Length must be a power of two");
            }
            this.m = (int)Math.Log(fftLength, 2.0);
            this.fftLength = fftLength;
            this.fftBuffer = new Complex[fftLength];
            this.fftArgs = new FftEventArgs(fftBuffer);
        }

        bool IsPowerOfTwo(int x)
        {
            return (x & (x - 1)) == 0;
        }


        public void Reset()
        {
            count = 0;
            maxValue = minValue = 0;
        }

        public void Add(float value)
        {
            if (PerformFFT)
            {
                fftBuffer[fftPos].X = (float)(value * FastFourierTransform.HammingWindow(fftPos,
fftBuffer.Length));
                fftBuffer[fftPos].Y = 0;
                fftPos++;
                if (fftPos >= fftBuffer.Length)
                {
                    fftPos = 0;
                    // 1024 = 2^10
                    FastFourierTransform.FFT(true, m, fftBuffer);

                    if(FftCalculated != null) FftCalculated(this, fftArgs);
                }
            }

            maxValue = Math.Max(maxValue, value);
            minValue = Math.Min(minValue, value);
            count++;
            if (count >= NotificationCount && NotificationCount > 0)
            {
                if (MaximumCalculated != null)
                {
                    MaximumCalculated(this, new MaxSampleEventArgs(minValue, maxValue));
                }
                Reset();
            }
        }
    }

    public class MaxSampleEventArgs : EventArgs
    {
        [DebuggerStepThrough]
        public MaxSampleEventArgs(float minValue, float maxValue)
        {
            this.MaxSample = maxValue;
            this.MinSample = minValue;
        }
        public float MaxSample { get; private set; }
        public float MinSample { get; private set; }
    }

    public class FftEventArgs : EventArgs
    {
        [DebuggerStepThrough]
        public FftEventArgs(Complex[] result)
        {
            this.Result = result;
        }
        public Complex[] Result { get; private set; }
    }
}


Тогда для определения частоты порции из первых 1024 сэмплов Wav-файла можно использовать
вот такой код:

using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Linq;
using System.Text;
using System.Windows.Forms;

using NAudio;
using NAudio.Wave;
using NAudio.Wave.SampleProviders;

namespace WindowsFormsTest1
{
    public partial class Form1 : Form
    {        

        private float parabolic(float[] f, int peak)
        {
            if (peak == 0) return f[0];
            var xv = 0.5f * (f[peak - 1] - f[peak + 1]) / (f[peak - 1] - 2 * f[peak]
+ f[peak + 1]) + peak;

            return xv;
        }       


        public Form1()
        {
            InitializeComponent();            
        }

        void PrintFrequency(float[] samples, int n_samples, WaveFormat fmt)
        {
            textBox1.Text = "";
            for (int i = 0; i < fmt.Channels; i++)
            {
                SampleAggregator aggregator = new SampleAggregator(n_samples);
                aggregator.PerformFFT = true;

                int j;
                float f;

                //выделяем данные одного канала
                for (j = 0; j < n_samples; j++)
                {
                    int index = (j * fmt.Channels) + i;
                    f = samples[index];
                    aggregator.Add(f);
                }

                float[] fft_x = new float[aggregator.FftBuffer.Length / 2]; //только
первая половина БПФ имеет смысл

                for (j = 0; j < fft_x.Length; j++)
                {
                   float real = aggregator.FftBuffer[j].X;
                   float imag = aggregator.FftBuffer[j].Y;
                   fft_x[j] = (float)Math.Sqrt(real * real + imag * imag); //получаем
амплитуду  
                }

                fft_x[0] = 0.0f;//избавляемся от постоянной составляющей

                int i_peak = fft_x.ToList().IndexOf(fft_x.Max());                

                for (j = 0; j < fft_x.Length; j++)
                {
                    fft_x[j] = (float)Math.Log(Math.Abs(aggregator.FftBuffer[j].X));
                }

                var i_interp = parabolic(fft_x, i_peak);

                float freq = fmt.SampleRate * i_interp / (float)n_samples;

                textBox1.Text += ("Channel " + i.ToString() + ": " + freq.ToString()
+ " Hz" + Environment.NewLine);
            }

        }


        private void button1_Click(object sender, EventArgs e)
        {

            WaveStream readerStream = new WaveFileReader("c:\\Test\\sound_01.wav");
            WaveStream pcmStream;
            WaveStream stream;

            //создаем поток в PCM-формате
            if (readerStream.WaveFormat.Encoding != WaveFormatEncoding.Pcm)
            {
                pcmStream = WaveFormatConversionStream.CreatePcmStream(readerStream);
                stream = new BlockAlignReductionStream(pcmStream);
            }
            else
            {
                pcmStream = readerStream;
                stream = readerStream;
            }

            float[] samples;
            const int N_SAMPLES = 1024; //количество сэмплов для спектрального анализа
            ISampleProvider prov;

            using(stream)
            using(readerStream)
            using (pcmStream)
            {
                prov = stream.ToSampleProvider();                
                samples = new float[N_SAMPLES * prov.WaveFormat.Channels];

                int res = prov.Read(samples, 0, N_SAMPLES * prov.WaveFormat.Channels);
                if (res < N_SAMPLES * prov.WaveFormat.Channels) throw new Exception("Not
enough data");

            }

            PrintFrequency(samples,N_SAMPLES,prov.WaveFormat);

        }  

    }

}

Какие есть алгоритмы шифрования изображений, устойчивые к трансформации?

#фурье #fft #изображения #wavelet #шифрование

Какие есть алгоритмы шифрования изображений, устойчивые к изменению размера изображения,
поворотам на углы, кратные 90 градусам?
Т.е. зашифровал исходное изображение - оно на себя не похоже, набор узоров, загрузил
на фотохостинг, где его в т.ч. уменьшили. Взял уменьшенный вариант, и тем же ключом
расшифровал в уменьшенную копию исходной картинки. Есть такие решения?
В идеале — вообще зашифровали - распечатали - сфотографировали на мобильник распечатку
- дешифровали в исходное, пусть такого же уменьшенного разрешения и плюс искажения
перспективы. Т.е. решение, не привязанное жестко к пикселам, и к точным значениям цветов/яркостей.
Upd. как шаг к анонимности интернетов в эпоху тотального контроля государств и корпораций,
имеет смысл добиваться такого вида зашифрованного изображения, что автоматическими
методами сложно провести его связь с исходником. Т.е. например, гистограмма зашифрованного
никак не должна кореллировать с гистограммой исходника. Например, если порезать картинку
на блоки и переставить их местами, гистограмма никак не изменится. Если канал инвертнуть,
гистограмма его просто отзеркалится. Поэтому, повторю, интересно услышать тех, кто
работал с FFT (быстрое преобразование Фурье) и вейвлетами.    

Ответы

Ответ 1

UPD [2013-02-07] Новый алгоритм «Chan»
Не оставляет меня в покое эта тема. Написал новый алгоритм шифрования.
Алгоритм работает очень долго, посему просьба не тестировать на больших картинках
и набраться терпения: кнопки «Encrypt» и «Decrypt» понимают с первого раза, не стреляйте
в них очередью.
Примеры:

Оригинал
Зашифрованное изображение
Расшифрованное изображение
Расшифрованное с неверным паролем

Основан на преобразовании цветов пикселей, но работает весьма неплохо, на мой взгляд. 
Недостатки:

шифрование слишком слабое; перебо́ров за 10–20 можно подобрать такой пароль, который
даст представление о содержании изображения (хотя цвета могут быть искажены)
JPEG-компрессия сказывается на расшифрованном изображении (появляется рябь), поэтому
на тестовой странице я выдаю результат в PNG
очень много коллизий ключей
очень медленно, необходимо оптимизировать

Положительные черты:

устойчивость к любым трансформациям (кроме цвета), будь то поворот/разворот или даже
посторонние данные, например:
Исходное: зашифрованная ранее картинка находится на нормальной (к примеру, фото плаката
с изображением зашифрованной фотографии)
Расшифрованное: зашифрованная картинка «распаковалась»

Есть идея добавить свойство (bool) loseless, в зависимости от которого алгоритм и,
соответственно, результат, будет меняться:

при loseless = on – будет работать как сейчас;
при loseless = off –  шифрование будет значительно надежнее, но с потерей двух цветов:
0xFFFFFF станет 0xFFFFFE, а 0x000000 – 0x000001 (то есть, и в зашифрованном и в дешифрованном
абс. белый и абс. черный будут отсутствовать).

Или можно не заморачиваться и прописать loseless = off как единственный вариант.
@sergiks, а у тебя как дела с этой задачей?
-----------------------------------------------------------------------------------------
UPD [2013-01-30] Новая версия алгоритма «мозаика»

Доделал вот прототип: http://image.lotoflot.com/mosaic_crypt.php
Выдерживает повороты кратные 90 град., но с изменением размера беда: очень сильные
искажения после дешифровки.
Удачный пример

Оригинал
Зашифрованное изображение
Расшифрованное изображение

Лавина хорошая благодаря md5. Ключ внутри функции раскладывается на составляющие,
каждый из которых хешируется. Для нормальной работы нужен массивный набор данных (а
не 2-3 символа), поэтому решил не изобретать и подключить md5.
Параметр Pieces определяет, на сколько кусков в обоих направлениях (X, Y) будет порезана
картинка. Чем больше значение, тем больше блоков. Обычно лучше устанавливать значение
в диапазоне 10—30 для картинок размера 300—2000 px.
Я как-то давно на PHP начал писать класс-прослойку для GD. Так вот эти функции шифровки/дешифровки
я запихнул в класс расширения этой прослойки, так что код будет слегка запутанным.
Но если все же заинтересует, могу поделиться.


Ответ 2

Конкретных алгоритмов не скажу, но наверняка они есть. Наверное надо копать в сторону
именно графической шифровки, к примеру поставить точки которые организуют квадрат чтобы
уйти от зависимости от ориентации картинки, и внутри этого генерировать данные которые
представляют исходную картинку. самое тупое - например точка 2x2 или 3x3 - 1, нет такой
точки - 0. Можно наверное более сложно придумать - яркая точка начала, от неё рисовать
примитивы - линия под углом - определённый набор битов (например угол задаётся исходным
куском битов, цвет каждой следующей линии - следующая пачка данных) ну и т.д. Ваша
задача - графически представить нужные вам данные - как - простор для фантазии. 
Так же можно использовать тупо набор известных вам однозначно обратимых фильтрах
(сдвиги, повороты и т.д.) которые Вы применяете в определённом порядке, а для расшифровки
- в обратном (например картинку сдвигает на 10 пикселов, поворачиваете на 15 градусов,
сдвигаете на 3 пиксела вниз, инвертируете цвета, сдвигаете цвета на 20, меняете каналы
местами и т.д.).


Ответ 3

Я бы посоветовал копать в сторону стеганографических трюков с изображениями, когда
модифицируется наименее значимый бит (LSB) пикселя или что-то подобное. Обзор некоторых
из таких методов здесь
Наиболее перспективный для вас метод на мой взгляд это redundant pattern encoding,
то есть когда шифровка размещается не по всему изображению, а дублируется в нескольких
частях (pattern), так что при обрезке изображения останутся кодированные куски из которых
можно будет восстановить сообщение.
Теперь рассмотрим способ восстановления исходного шифросообщения после изменения
размера изображения (паттерна). Здесь может помочь один из алгоритмов избыточного кодирования,
например метод Хамминга или что-то подобное. Благо таких способов довольно много, надо
просто задаться процентом сколько бит надо уметь восстанавливать из потерянных - это
чуть ли не целая отрасль математики.


Ответ 4

Почему вы так быстро отказались от шифрование данных каждого пикселя независимо друг
от друга? Изменения размеров и артефакты сжатия такому методу не страшны. Я, обдумывая
ту же идею, склонен предполагать, что можно найти достаточно устойчивый (к взлому)
шифр. Но так как в области криптографии моих знаний ноль, с этим вопросом я и пришёл сюда.

Нормализация кросс-корреляции

Доброго времени суток! Я пытаюсь написать программу для поиска шаблона в сигнале. Сигнал переодический, квазистационарный. Задача - получить адекватный коэффициент корреляции(КК). Перед основными вычислениями произвожу удаление среднего из выборок: x = x - mean(x) Для расчета КК использую подход на основе БПФ: corr(x, y) = F'( F(x) * conj(F(y)) ), где F() - прямое БПФ, F'() - обратное БПФ, conj() - получение комплексно сопряженного. Для того чтобы работало БПФ использую дополнение нулями исходных массивов x и y до степени 2 следующим образом [000000xxxxxxx000] [yyyyyyy000000000] После выполнения процедуры corr() получаю кросс-корреляционную функцию, максимум которой вроде как является искомым КК, а позиция максимума соответствует сдвигу шаблона в сигнале. Однако, чтобы получить КК полученный максимум надо нормировать(чтобы получить значение в диапазоне от -1 до 1) Вопрос собственно в том как правильно нормировать полученный максимум? P.S. если будет необходимо могу привести исходный код.

Ответ

Хм. IronVbif, а ведь Вы похоже правы. Система дискретна.Можно либо применить критерий Пирсона к начальной задаче (и сразу получить результат). А также можно получить нормированную корреляцию путем деления полученной функции на корень квадратный из дисперсии: normcov=f(x)/(sqrt((sum((f(x) - mean(x))^2)/n)). Похоже, это и есть ответ.

Функция вызывается с неправильными значениями переменных

Суть вопроса такова. Есть такой код:
#include "stdafx.h"
#define PI 3.1415926535897932384626433832795 #define FR 3 #define SAMPLES 1024
int _tmain(int argc, _TCHAR* argv[]) { base signal[SAMPLES];
for (int i = 0; i < SAMPLES; ++i) { signal[i] = (sin(2*PI*i*FR), 0); }
FFTCalculate(signal, SAMPLES, false);
FILE* f_pointer = fopen("fft_result.txt", "w");
for (int i = 0; i < SAMPLES; ++i) { char num = (char)signal[i].real(); int out = fprintf(f_pointer, &num); } int out = fclose(f_pointer); return 0; }
в нем вызывается функция FFTCalculate(signal, SAMPLES, false). Ниже приведен код функции:
void FFTCalculate(base signal[], int n, bool invert) { int log_N; double x = frexp((double)n, &log_N);
calc_rev(n, log_N);
for(int i = 0; i < n; ++i) if (i < rev[i]) swap(signal[i], signal[rev[i]]);
for(int len = 2; len <= n; len<<=1) { double ang = 2*PI/len * (invert?-1:1); int len2 = len>>1;
base wlen (cos(ang), sin(ang)); wlen_pw[0] = base(1, 0); for (int i=0; i for (int i=0; i for(; pu != pu_end; ++pu, ++pv, ++pu){ t = *pv * *pw; *pv = *pu - t; *pu += t; } } } if(invert) for (int i=0; iПри запуске программы появляется ошибка, связанная с выходом за пределы массива. Начал дебажить и выяснил, что когда вызывается функция FFTCalculate(signal, SAMPLES, false), то ей почему то передаются не signal, SAMPLES = 1024 и false, а sugnal, 1245452 и true см. скриншот ниже. В чем дело?

Ответ

Я проверил на GCC 4.8.2 под linux. Параметры FFTCalculate передаются корректно. Думаю что, то поведение которое ты видишь в Visual Studio является unexpected behavior. Твоя программа под линуском сразу упала с seg fault.
Ниже пример кода с ошибкой. Ты не проверяешь значения rev[i], что оно может быть больше n (1024).
calc_rev(n, log_N);
for(int i = 0; i < n; ++i) if (i < rev[i]) swap(signal[i], signal[rev[i]]);
Вот что у меня вывелось для первых десяти итераций из массива rev: 0 1024 1536 1280 1792 1152 1664 1408 1920
Надеюсь помог.

Алгоритмы активного шумоподавления

Меня интересуют алгоритмы активного шумоподавления. Углубляясь в детали реализации, выполняется спектральное вычитание "на ходу".
Имеется некоторый профиль шума и небольшой длины запись, из которой этот шум нужно удалить ("вычесть").
Так у меня вопрос по конкретному алгоритму спектрального вычитания: как оно выполняется и, собственно, над чем? В общих чертах я себе это представляю как вычитание из спектра (соответствие частот амплитудам гармоник) сигнала спектра шума и обратное преобразование Фурье. Вообще, вроде как, применяемое здесь FFT даёт не соответствие частот амплитудам, а просто массив комплексных чисел. Как тут вычитать что из чего не представляю...
P.S.: Не нужно давать ссылки на курсы по обработке звука и спектральному анализу сигналов. Просто вопрос - просто ясный ответ.

Есть такие специализированные программы, как, например, Noise Gator. Так вот, если кому известно, подскажите, какие алгоритмы шумоподавления ими используются (в частности, используется ли спектральное вычитание; если да, то откуда программа берёт слепок шума).

Ответ

Алгоритмы удаления шума невозможно рассмотреть без углубления в математику, акустику и теорию спектрального анализа сигналов.
Любой анализ сигнала с помощью преобразований Фурье предполагает, что сигнал является стационарным на исследуемом отрезке. Поэтому, если сигнал нестационарный, он разбивается на отдельные отрезки, называемые окнами. Выбор размера окна зависит от типа исследуемого сигнала, обычно около 25 - 50 мс для звукового сигнала (меньшие значения - для человеческой речи, большие - для музыки, особенно состоящей из струнных смычковых инструментов). Можно использовать перекрывающиеся окна, для повышения точности анализа.
Однако, просто так применить преобразование Фурье к обрезанным окнам нельзя, при этом некорректно обрабатываются граничные области отрезков. Для решения этой проблемы сигнал предварительно умножают на специальную весовую функцию ("оконную"). Примеры оконных функций см. в статье Оконное преобразование Фурье
Далее выполняется непосредственно преобразование Фурье. Оно дает в результате спектр сигнала, т.е. значение комплексной амплитуды для различных диапазонов частот. Из него и надо вычитать спектр шума. Из модуля комплексной амплитуды сигнала вычитается модуль комплексной амплитуды шума, умноженный на некий коэффициент; если результат отрицательный, он заменяется на ноль. Фазовый компонент оставляется нетронутым. К результату вычитания можно применить обратное преобразование Фурье, и получить "очищенный" сигнал.
Итоговый алгоритм шумоподавления при наличии известного образца шума:
Разделение сигнала на окна Применение оконного преобразования Фурье к окнам Вычитание (по модулю) спектра амплитуды шума из спектра амплитуды сигнала:
A = Max ( A с. - k * А ш. ; 0)
где k - коэффициент, подбираемый опытным путем
Применение обратного преобразования Фурье к результату
Размер окна, перекрытие окон, тип применяемой оконной функции подбираются опытным путем.
Ссылки
Removing noise from audio using Fourier transform in Matlab
How can I select an optimal window for Short Time Fourier Transform?
How to select frequency resolution and window size in FFT?