Merge branch 'spectrogramoptim' of https://github.com/jdpurcell/subtitleedit into playselection

This commit is contained in:
J.D. Purcell 2015-09-02 23:41:32 -04:00
commit 2c365d078f
5 changed files with 956 additions and 374 deletions

View File

@ -8,13 +8,21 @@ namespace Nikse.SubtitleEdit.Core
{
unsafe public class FastBitmap
{
private struct PixelData
public struct PixelData
{
public byte Blue;
public byte Green;
public byte Red;
public byte Alpha;
public PixelData(Color c)
{
Alpha = c.A;
Red = c.R;
Green = c.G;
Blue = c.B;
}
public override string ToString()
{
return "(" + Alpha + ", " + Red + ", " + Green + ", " + Blue + ")";
@ -82,6 +90,12 @@ namespace Nikse.SubtitleEdit.Core
data->Blue = color.B;
}
public void SetPixel(int x, int y, PixelData color)
{
var data = (PixelData*)(_pBase + y * _width + x * sizeof(PixelData));
*data = color;
}
public void SetPixel(int x, int y, Color color, int length)
{
var data = (PixelData*)(_pBase + y * _width + x * sizeof(PixelData));

View File

@ -1,162 +0,0 @@
using System;
namespace Nikse.SubtitleEdit.Core
{
/// <summary>
/// Fourier transform
///
/// *****************************************************************************
/// *
/// * Copyright (c) 2002, Wilhelm Kurz. All Rights Reserved.
/// * wkurz@foni.net
/// *
/// * This file is provided for demonstration and educational uses only.
/// * Permission to use, copy, modify and distribute this file for
/// * any purpose and without fee is hereby granted.
/// *
/// *****************************************************************************
/// Converted/optimized by Nikse from vb code: http://www.wilhelm-kurz-software.de/dynaplot/applicationnotes/spectrogram.htm
/// </summary>
internal class Fourier
{
public const double W0Hanning = 0.5;
public const double W0Hamming = 0.54;
public const double W0Blackman = 0.42;
private const double Pi = 3.14159265358979;
private double[] cosarray;
private double[] sinarray;
private bool _forward;
private int _arraySize;
private int _ldArraysize = 0;
public Fourier(int arraySize, bool forward)
{
_arraySize = arraySize;
_forward = forward;
cosarray = new double[arraySize];
sinarray = new double[arraySize];
double sign = 1.0;
if (forward)
sign = -1.0;
double phase0 = 2.0 * Pi / arraySize;
for (int i = 0; i <= arraySize - 1; i++)
{
sinarray[i] = sign * Math.Sin(phase0 * i);
cosarray[i] = Math.Cos(phase0 * i);
}
int j = _arraySize;
while (j != 1)
{
_ldArraysize++;
j /= 2;
}
}
public void MagnitudeSpectrum(double[] real, double[] imag, double w0, double[] magnitude)
{
int i;
magnitude[0] = Math.Sqrt(SquareSum(real[0], imag[0]));
for (i = 1; i <= (_arraySize / 2 - 1); i++)
magnitude[i] = (Math.Sqrt(SquareSum(real[i], imag[i]) + SquareSum(real[_arraySize - i], imag[_arraySize - i]))) / w0;
}
public static double Hanning(int n, int j)
{
return W0Hanning - 0.5 * Math.Cos(2.0 * Pi * j / n);
}
public static double Hamming(int n, int j)
{
return W0Hamming - 0.46 * Math.Cos(2.0 * Pi * j / n);
}
public static double Blackman(int n, int j)
{
return W0Blackman - 0.5 * Math.Cos(2.0 * Pi * j / n) + 0.08 * Math.Cos(4.0 * Pi * j / n);
}
private static void Swap(ref double a, ref double b)
{
double temp = a;
a = b;
b = temp;
}
private static double SquareSum(double a, double b)
{
return a * a + b * b;
}
public void FourierTransform(double[] real, double[] imag)
{
int i;
if (_forward)
{
for (i = 0; i <= _arraySize - 1; i++)
{
real[i] /= _arraySize;
imag[i] /= _arraySize;
}
}
int k;
int j = 0;
for (i = 0; i <= _arraySize - 2; i++)
{
if (i < j)
{
Swap(ref real[i], ref real[j]);
Swap(ref imag[i], ref imag[j]);
}
k = _arraySize / 2;
while (k <= j)
{
j -= k;
k /= 2;
}
j += k;
}
int a = 2;
int b = 1;
for (int count = 1; count <= _ldArraysize; count++)
{
int c0 = _arraySize / a;
int c1 = 0;
for (k = 0; k <= b - 1; k++)
{
i = k;
while (i < _arraySize)
{
int arg = i + b;
double prodreal;
double prodimag;
if (k == 0)
{
prodreal = real[arg];
prodimag = imag[arg];
}
else
{
prodreal = real[arg] * cosarray[c1] - imag[arg] * sinarray[c1];
prodimag = real[arg] * sinarray[c1] + imag[arg] * cosarray[c1];
}
real[arg] = real[i] - prodreal;
imag[arg] = imag[i] - prodimag;
real[i] += prodreal;
imag[i] += prodimag;
i += a;
}
c1 += c0;
}
a *= 2;
b *= 2;
}
}
}
}

View File

@ -132,7 +132,6 @@
<Compile Include="Forms\RemoveTextForHI.cs" />
<Compile Include="Forms\RemoveTextForHISettings.cs" />
<Compile Include="Forms\SplitLongLinesHelper.cs" />
<Compile Include="Fourier.cs" />
<Compile Include="HistoryItem.cs" />
<Compile Include="HtmlUtil.cs" />
<Compile Include="IfoParser.cs" />
@ -148,6 +147,7 @@
<Compile Include="NoBreakAfterItem.cs" />
<Compile Include="Paragraph.cs" />
<Compile Include="Properties\AssemblyInfo.cs" />
<Compile Include="RealFFT.cs" />
<Compile Include="RegistryUtil.cs" />
<Compile Include="Settings.cs" />
<Compile Include="SpellCheckWord.cs" />

666
libse/RealFFT.cs Normal file
View File

@ -0,0 +1,666 @@
// Copyright(C) 1996-2001 Takuya Ooura
// C# port by J.D. Purcell
// You may use, copy, modify this code for any purpose and without fee.
using System;
namespace Nikse.SubtitleEdit.Core
{
public class RealFFT
{
private int _length;
private int[] _ip;
private double[] _w;
public readonly double ForwardScaleFactor;
public readonly double ReverseScaleFactor;
public RealFFT(int length)
{
if (length < 2 || (length & (length - 1)) != 0)
{
throw new ArgumentException("length", "FFT length must be at least 2 and a power of 2.");
}
_length = length;
_ip = new int[2 + (1 << (Convert.ToInt32(Math.Log(Math.Max(length / 4, 1), 2)) / 2))];
_w = new double[length / 2];
ForwardScaleFactor = length;
ReverseScaleFactor = 0.5;
}
public void ComputeForward(double[] buff)
{
Compute(buff, false);
}
public void ComputeReverse(double[] buff)
{
Compute(buff, true);
}
private void Compute(double[] buff, bool reverse)
{
if (buff.Length < _length)
{
throw new ArgumentException("buff", "Buffer length must be greater than or equal to the FFT length.");
}
rdft(_length, reverse, buff, _ip, _w);
}
private static void rdft(int n, bool rev, double[] a, int[] ip, double[] w)
{
int nw, nc;
double xi;
nw = ip[0];
if (n > (nw << 2))
{
nw = n >> 2;
makewt(nw, ip, w);
}
nc = ip[1];
if (n > (nc << 2))
{
nc = n >> 2;
makect(nc, ip, w, nw);
}
if (!rev)
{
if (n > 4)
{
bitrv2(n, ip, a);
cftfsub(n, a, w);
rftfsub(n, a, nc, w, nw);
}
else if (n == 4)
{
cftfsub(n, a, w);
}
xi = a[0] - a[1];
a[0] += a[1];
a[1] = xi;
}
else
{
a[1] = 0.5 * (a[0] - a[1]);
a[0] -= a[1];
if (n > 4)
{
rftbsub(n, a, nc, w, nw);
bitrv2(n, ip, a);
cftbsub(n, a, w);
}
else if (n == 4)
{
cftfsub(n, a, w);
}
}
}
/* -------- initializing routines -------- */
private static void makewt(int nw, int[] ip, double[] w)
{
int j, nwh;
double delta, x, y;
ip[0] = nw;
ip[1] = 1;
if (nw > 2)
{
nwh = nw >> 1;
delta = Math.Atan(1.0) / nwh;
w[0] = 1;
w[1] = 0;
w[nwh] = Math.Cos(delta * nwh);
w[nwh + 1] = w[nwh];
if (nwh > 2)
{
for (j = 2; j < nwh; j += 2)
{
x = Math.Cos(delta * j);
y = Math.Sin(delta * j);
w[j] = x;
w[j + 1] = y;
w[nw - j] = y;
w[nw - j + 1] = x;
}
bitrv2(nw, ip, w);
}
}
}
private static void makect(int nc, int[] ip, double[] c, int nw)
{
int j, nch;
double delta;
ip[1] = nc;
if (nc > 1)
{
nch = nc >> 1;
delta = Math.Atan(1.0) / nch;
c[nw] = Math.Cos(delta * nch);
c[nw + nch] = 0.5 * c[nw];
for (j = 1; j < nch; j++)
{
c[nw + j] = 0.5 * Math.Cos(delta * j);
c[nw + nc - j] = 0.5 * Math.Sin(delta * j);
}
}
}
/* -------- child routines -------- */
private static void bitrv2(int n, int[] ip, double[] a)
{
int j, j1, k, k1, l, m, m2;
double xr, xi, yr, yi;
ip[2] = 0;
l = n;
m = 1;
while ((m << 3) < l)
{
l >>= 1;
for (j = 0; j < m; j++)
{
ip[m + j + 2] = ip[j + 2] + l;
}
m <<= 1;
}
m2 = 2 * m;
if ((m << 3) == l)
{
for (k = 0; k < m; k++)
{
for (j = 0; j < k; j++)
{
j1 = 2 * j + ip[k + 2];
k1 = 2 * k + ip[j + 2];
xr = a[j1];
xi = a[j1 + 1];
yr = a[k1];
yi = a[k1 + 1];
a[j1] = yr;
a[j1 + 1] = yi;
a[k1] = xr;
a[k1 + 1] = xi;
j1 += m2;
k1 += 2 * m2;
xr = a[j1];
xi = a[j1 + 1];
yr = a[k1];
yi = a[k1 + 1];
a[j1] = yr;
a[j1 + 1] = yi;
a[k1] = xr;
a[k1 + 1] = xi;
j1 += m2;
k1 -= m2;
xr = a[j1];
xi = a[j1 + 1];
yr = a[k1];
yi = a[k1 + 1];
a[j1] = yr;
a[j1 + 1] = yi;
a[k1] = xr;
a[k1 + 1] = xi;
j1 += m2;
k1 += 2 * m2;
xr = a[j1];
xi = a[j1 + 1];
yr = a[k1];
yi = a[k1 + 1];
a[j1] = yr;
a[j1 + 1] = yi;
a[k1] = xr;
a[k1 + 1] = xi;
}
j1 = 2 * k + m2 + ip[k + 2];
k1 = j1 + m2;
xr = a[j1];
xi = a[j1 + 1];
yr = a[k1];
yi = a[k1 + 1];
a[j1] = yr;
a[j1 + 1] = yi;
a[k1] = xr;
a[k1 + 1] = xi;
}
}
else
{
for (k = 1; k < m; k++)
{
for (j = 0; j < k; j++)
{
j1 = 2 * j + ip[k + 2];
k1 = 2 * k + ip[j + 2];
xr = a[j1];
xi = a[j1 + 1];
yr = a[k1];
yi = a[k1 + 1];
a[j1] = yr;
a[j1 + 1] = yi;
a[k1] = xr;
a[k1 + 1] = xi;
j1 += m2;
k1 += m2;
xr = a[j1];
xi = a[j1 + 1];
yr = a[k1];
yi = a[k1 + 1];
a[j1] = yr;
a[j1 + 1] = yi;
a[k1] = xr;
a[k1 + 1] = xi;
}
}
}
}
private static void cftfsub(int n, double[] a, double[] w)
{
int j, j1, j2, j3, l;
double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
l = 2;
if (n > 8)
{
cft1st(n, a, w);
l = 8;
while ((l << 2) < n)
{
cftmdl(n, l, a, w);
l <<= 2;
}
}
if ((l << 2) == n)
{
for (j = 0; j < l; j += 2)
{
j1 = j + l;
j2 = j1 + l;
j3 = j2 + l;
x0r = a[j] + a[j1];
x0i = a[j + 1] + a[j1 + 1];
x1r = a[j] - a[j1];
x1i = a[j + 1] - a[j1 + 1];
x2r = a[j2] + a[j3];
x2i = a[j2 + 1] + a[j3 + 1];
x3r = a[j2] - a[j3];
x3i = a[j2 + 1] - a[j3 + 1];
a[j] = x0r + x2r;
a[j + 1] = x0i + x2i;
a[j2] = x0r - x2r;
a[j2 + 1] = x0i - x2i;
a[j1] = x1r - x3i;
a[j1 + 1] = x1i + x3r;
a[j3] = x1r + x3i;
a[j3 + 1] = x1i - x3r;
}
}
else
{
for (j = 0; j < l; j += 2)
{
j1 = j + l;
x0r = a[j] - a[j1];
x0i = a[j + 1] - a[j1 + 1];
a[j] += a[j1];
a[j + 1] += a[j1 + 1];
a[j1] = x0r;
a[j1 + 1] = x0i;
}
}
}
private static void cftbsub(int n, double[] a, double[] w)
{
int j, j1, j2, j3, l;
double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
l = 2;
if (n > 8)
{
cft1st(n, a, w);
l = 8;
while ((l << 2) < n)
{
cftmdl(n, l, a, w);
l <<= 2;
}
}
if ((l << 2) == n)
{
for (j = 0; j < l; j += 2)
{
j1 = j + l;
j2 = j1 + l;
j3 = j2 + l;
x0r = a[j] + a[j1];
x0i = -a[j + 1] - a[j1 + 1];
x1r = a[j] - a[j1];
x1i = -a[j + 1] + a[j1 + 1];
x2r = a[j2] + a[j3];
x2i = a[j2 + 1] + a[j3 + 1];
x3r = a[j2] - a[j3];
x3i = a[j2 + 1] - a[j3 + 1];
a[j] = x0r + x2r;
a[j + 1] = x0i - x2i;
a[j2] = x0r - x2r;
a[j2 + 1] = x0i + x2i;
a[j1] = x1r - x3i;
a[j1 + 1] = x1i - x3r;
a[j3] = x1r + x3i;
a[j3 + 1] = x1i + x3r;
}
}
else
{
for (j = 0; j < l; j += 2)
{
j1 = j + l;
x0r = a[j] - a[j1];
x0i = -a[j + 1] + a[j1 + 1];
a[j] += a[j1];
a[j + 1] = -a[j + 1] - a[j1 + 1];
a[j1] = x0r;
a[j1 + 1] = x0i;
}
}
}
private static void cft1st(int n, double[] a, double[] w)
{
int j, k1, k2;
double wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
x0r = a[0] + a[2];
x0i = a[1] + a[3];
x1r = a[0] - a[2];
x1i = a[1] - a[3];
x2r = a[4] + a[6];
x2i = a[5] + a[7];
x3r = a[4] - a[6];
x3i = a[5] - a[7];
a[0] = x0r + x2r;
a[1] = x0i + x2i;
a[4] = x0r - x2r;
a[5] = x0i - x2i;
a[2] = x1r - x3i;
a[3] = x1i + x3r;
a[6] = x1r + x3i;
a[7] = x1i - x3r;
wk1r = w[2];
x0r = a[8] + a[10];
x0i = a[9] + a[11];
x1r = a[8] - a[10];
x1i = a[9] - a[11];
x2r = a[12] + a[14];
x2i = a[13] + a[15];
x3r = a[12] - a[14];
x3i = a[13] - a[15];
a[8] = x0r + x2r;
a[9] = x0i + x2i;
a[12] = x2i - x0i;
a[13] = x0r - x2r;
x0r = x1r - x3i;
x0i = x1i + x3r;
a[10] = wk1r * (x0r - x0i);
a[11] = wk1r * (x0r + x0i);
x0r = x3i + x1r;
x0i = x3r - x1i;
a[14] = wk1r * (x0i - x0r);
a[15] = wk1r * (x0i + x0r);
k1 = 0;
for (j = 16; j < n; j += 16)
{
k1 += 2;
k2 = 2 * k1;
wk2r = w[k1];
wk2i = w[k1 + 1];
wk1r = w[k2];
wk1i = w[k2 + 1];
wk3r = wk1r - 2 * wk2i * wk1i;
wk3i = 2 * wk2i * wk1r - wk1i;
x0r = a[j] + a[j + 2];
x0i = a[j + 1] + a[j + 3];
x1r = a[j] - a[j + 2];
x1i = a[j + 1] - a[j + 3];
x2r = a[j + 4] + a[j + 6];
x2i = a[j + 5] + a[j + 7];
x3r = a[j + 4] - a[j + 6];
x3i = a[j + 5] - a[j + 7];
a[j] = x0r + x2r;
a[j + 1] = x0i + x2i;
x0r -= x2r;
x0i -= x2i;
a[j + 4] = wk2r * x0r - wk2i * x0i;
a[j + 5] = wk2r * x0i + wk2i * x0r;
x0r = x1r - x3i;
x0i = x1i + x3r;
a[j + 2] = wk1r * x0r - wk1i * x0i;
a[j + 3] = wk1r * x0i + wk1i * x0r;
x0r = x1r + x3i;
x0i = x1i - x3r;
a[j + 6] = wk3r * x0r - wk3i * x0i;
a[j + 7] = wk3r * x0i + wk3i * x0r;
wk1r = w[k2 + 2];
wk1i = w[k2 + 3];
wk3r = wk1r - 2 * wk2r * wk1i;
wk3i = 2 * wk2r * wk1r - wk1i;
x0r = a[j + 8] + a[j + 10];
x0i = a[j + 9] + a[j + 11];
x1r = a[j + 8] - a[j + 10];
x1i = a[j + 9] - a[j + 11];
x2r = a[j + 12] + a[j + 14];
x2i = a[j + 13] + a[j + 15];
x3r = a[j + 12] - a[j + 14];
x3i = a[j + 13] - a[j + 15];
a[j + 8] = x0r + x2r;
a[j + 9] = x0i + x2i;
x0r -= x2r;
x0i -= x2i;
a[j + 12] = -wk2i * x0r - wk2r * x0i;
a[j + 13] = -wk2i * x0i + wk2r * x0r;
x0r = x1r - x3i;
x0i = x1i + x3r;
a[j + 10] = wk1r * x0r - wk1i * x0i;
a[j + 11] = wk1r * x0i + wk1i * x0r;
x0r = x1r + x3i;
x0i = x1i - x3r;
a[j + 14] = wk3r * x0r - wk3i * x0i;
a[j + 15] = wk3r * x0i + wk3i * x0r;
}
}
private static void cftmdl(int n, int l, double[] a, double[] w)
{
int j, j1, j2, j3, k, k1, k2, m, m2;
double wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
m = l << 2;
for (j = 0; j < l; j += 2)
{
j1 = j + l;
j2 = j1 + l;
j3 = j2 + l;
x0r = a[j] + a[j1];
x0i = a[j + 1] + a[j1 + 1];
x1r = a[j] - a[j1];
x1i = a[j + 1] - a[j1 + 1];
x2r = a[j2] + a[j3];
x2i = a[j2 + 1] + a[j3 + 1];
x3r = a[j2] - a[j3];
x3i = a[j2 + 1] - a[j3 + 1];
a[j] = x0r + x2r;
a[j + 1] = x0i + x2i;
a[j2] = x0r - x2r;
a[j2 + 1] = x0i - x2i;
a[j1] = x1r - x3i;
a[j1 + 1] = x1i + x3r;
a[j3] = x1r + x3i;
a[j3 + 1] = x1i - x3r;
}
wk1r = w[2];
for (j = m; j < l + m; j += 2)
{
j1 = j + l;
j2 = j1 + l;
j3 = j2 + l;
x0r = a[j] + a[j1];
x0i = a[j + 1] + a[j1 + 1];
x1r = a[j] - a[j1];
x1i = a[j + 1] - a[j1 + 1];
x2r = a[j2] + a[j3];
x2i = a[j2 + 1] + a[j3 + 1];
x3r = a[j2] - a[j3];
x3i = a[j2 + 1] - a[j3 + 1];
a[j] = x0r + x2r;
a[j + 1] = x0i + x2i;
a[j2] = x2i - x0i;
a[j2 + 1] = x0r - x2r;
x0r = x1r - x3i;
x0i = x1i + x3r;
a[j1] = wk1r * (x0r - x0i);
a[j1 + 1] = wk1r * (x0r + x0i);
x0r = x3i + x1r;
x0i = x3r - x1i;
a[j3] = wk1r * (x0i - x0r);
a[j3 + 1] = wk1r * (x0i + x0r);
}
k1 = 0;
m2 = 2 * m;
for (k = m2; k < n; k += m2)
{
k1 += 2;
k2 = 2 * k1;
wk2r = w[k1];
wk2i = w[k1 + 1];
wk1r = w[k2];
wk1i = w[k2 + 1];
wk3r = wk1r - 2 * wk2i * wk1i;
wk3i = 2 * wk2i * wk1r - wk1i;
for (j = k; j < l + k; j += 2)
{
j1 = j + l;
j2 = j1 + l;
j3 = j2 + l;
x0r = a[j] + a[j1];
x0i = a[j + 1] + a[j1 + 1];
x1r = a[j] - a[j1];
x1i = a[j + 1] - a[j1 + 1];
x2r = a[j2] + a[j3];
x2i = a[j2 + 1] + a[j3 + 1];
x3r = a[j2] - a[j3];
x3i = a[j2 + 1] - a[j3 + 1];
a[j] = x0r + x2r;
a[j + 1] = x0i + x2i;
x0r -= x2r;
x0i -= x2i;
a[j2] = wk2r * x0r - wk2i * x0i;
a[j2 + 1] = wk2r * x0i + wk2i * x0r;
x0r = x1r - x3i;
x0i = x1i + x3r;
a[j1] = wk1r * x0r - wk1i * x0i;
a[j1 + 1] = wk1r * x0i + wk1i * x0r;
x0r = x1r + x3i;
x0i = x1i - x3r;
a[j3] = wk3r * x0r - wk3i * x0i;
a[j3 + 1] = wk3r * x0i + wk3i * x0r;
}
wk1r = w[k2 + 2];
wk1i = w[k2 + 3];
wk3r = wk1r - 2 * wk2r * wk1i;
wk3i = 2 * wk2r * wk1r - wk1i;
for (j = k + m; j < l + (k + m); j += 2)
{
j1 = j + l;
j2 = j1 + l;
j3 = j2 + l;
x0r = a[j] + a[j1];
x0i = a[j + 1] + a[j1 + 1];
x1r = a[j] - a[j1];
x1i = a[j + 1] - a[j1 + 1];
x2r = a[j2] + a[j3];
x2i = a[j2 + 1] + a[j3 + 1];
x3r = a[j2] - a[j3];
x3i = a[j2 + 1] - a[j3 + 1];
a[j] = x0r + x2r;
a[j + 1] = x0i + x2i;
x0r -= x2r;
x0i -= x2i;
a[j2] = -wk2i * x0r - wk2r * x0i;
a[j2 + 1] = -wk2i * x0i + wk2r * x0r;
x0r = x1r - x3i;
x0i = x1i + x3r;
a[j1] = wk1r * x0r - wk1i * x0i;
a[j1 + 1] = wk1r * x0i + wk1i * x0r;
x0r = x1r + x3i;
x0i = x1i - x3r;
a[j3] = wk3r * x0r - wk3i * x0i;
a[j3 + 1] = wk3r * x0i + wk3i * x0r;
}
}
}
private static void rftfsub(int n, double[] a, int nc, double[] c, int nw)
{
int j, k, kk, ks, m;
double wkr, wki, xr, xi, yr, yi;
m = n >> 1;
ks = 2 * nc / m;
kk = 0;
for (j = 2; j < m; j += 2)
{
k = n - j;
kk += ks;
wkr = 0.5 - c[nw + nc - kk];
wki = c[nw + kk];
xr = a[j] - a[k];
xi = a[j + 1] + a[k + 1];
yr = wkr * xr - wki * xi;
yi = wkr * xi + wki * xr;
a[j] -= yr;
a[j + 1] -= yi;
a[k] += yr;
a[k + 1] -= yi;
}
}
private static void rftbsub(int n, double[] a, int nc, double[] c, int nw)
{
int j, k, kk, ks, m;
double wkr, wki, xr, xi, yr, yi;
a[1] = -a[1];
m = n >> 1;
ks = 2 * nc / m;
kk = 0;
for (j = 2; j < m; j += 2)
{
k = n - j;
kk += ks;
wkr = 0.5 - c[nw + nc - kk];
wki = c[nw + kk];
xr = a[j] - a[k];
xi = a[j + 1] + a[k + 1];
yr = wkr * xr + wki * xi;
yi = wkr * xi - wki * xr;
a[j] -= yr;
a[j + 1] = yi - a[j + 1];
a[k] += yr;
a[k + 1] = yi - a[k + 1];
}
a[m + 1] = -a[m + 1];
}
}
}

View File

@ -124,6 +124,14 @@ namespace Nikse.SubtitleEdit.Core
}
}
public long LengthInSamples
{
get
{
return DataChunkSize / BlockAlign;
}
}
internal void WriteHeader(Stream toStream, int sampleRate, int numberOfChannels, int bitsPerSample, int dataSize)
{
const int fmtChunckSize = 16;
@ -219,7 +227,7 @@ namespace Nikse.SubtitleEdit.Core
{
PeaksPerSecond = peaksPerSecond;
ReadSampleDataValueDelegate readSampleDataValue = GetSampleDataRerader();
ReadSampleDataValueDelegate readSampleDataValue = GetSampleDataReader();
DataMinValue = int.MaxValue;
DataMaxValue = int.MinValue;
PeakSamples = new List<int>();
@ -258,7 +266,7 @@ namespace Nikse.SubtitleEdit.Core
public void GenerateAllSamples()
{
// determine how to read sample values
ReadSampleDataValueDelegate readSampleDataValue = GetSampleDataRerader();
ReadSampleDataValueDelegate readSampleDataValue = GetSampleDataReader();
// load data
_data = new byte[Header.DataChunkSize];
@ -325,7 +333,9 @@ namespace Nikse.SubtitleEdit.Core
private int ReadValue16Bit(ref int index)
{
int result = BitConverter.ToInt16(_data, index);
int result = (short)(
(_data[index ] ) |
(_data[index + 1] << 8));
index += 2;
return result;
}
@ -353,7 +363,7 @@ namespace Nikse.SubtitleEdit.Core
/// Determine how to read sample values
/// </summary>
/// <returns>Sample data reader that matches bits per sample</returns>
private ReadSampleDataValueDelegate GetSampleDataRerader()
private ReadSampleDataValueDelegate GetSampleDataReader()
{
ReadSampleDataValueDelegate readSampleDataValue;
switch (Header.BitsPerSample)
@ -392,115 +402,98 @@ namespace Nikse.SubtitleEdit.Core
public List<Bitmap> GenerateFourierData(int nfft, string spectrogramDirectory, int delayInMilliseconds)
{
const int bitmapWidth = 1024;
var bitmaps = new List<Bitmap>();
// setup fourier transformation
var f = new Fourier(nfft, true);
double divider = 2.0;
for (int k = 0; k < Header.BitsPerSample - 2; k++)
divider *= 2;
List<Bitmap> bitmaps = new List<Bitmap>();
SpectrogramDrawer drawer = new SpectrogramDrawer(nfft);
ReadSampleDataValueDelegate readSampleDataValue = GetSampleDataReader();
double sampleScale = 1.0 / (Math.Pow(2.0, Header.BitsPerSample - 1) * Header.NumberOfChannels);
// determine how to read sample values
ReadSampleDataValueDelegate readSampleDataValue = GetSampleDataRerader();
// set up one column of the spectrogram
var palette = new Color[nfft];
if (Configuration.Settings.VideoControls.SpectrogramAppearance == "Classic")
{
for (int colorIndex = 0; colorIndex < nfft; colorIndex++)
palette[colorIndex] = PaletteValue(colorIndex, nfft);
}
else
{
var list = SmoothColors(0, 0, 0, Configuration.Settings.VideoControls.WaveformColor.R,
Configuration.Settings.VideoControls.WaveformColor.G,
Configuration.Settings.VideoControls.WaveformColor.B, nfft);
for (int i = 0; i < nfft; i++)
palette[i] = list[i];
}
// read sample values
DataMinValue = int.MaxValue;
DataMaxValue = int.MinValue;
var samples = new List<int>();
int index = 0;
int sampleSize = nfft * bitmapWidth;
int count = 0;
long totalSamples = 0;
// write delay (if any)
int delaySampleCount = (int)(Header.SampleRate * (delayInMilliseconds / TimeCode.BaseUnit));
for (int i = 0; i < delaySampleCount; i++)
// other code (e.g. generating peaks) doesn't handle negative delays, so we'll do the same for now
delaySampleCount = Math.Max(delaySampleCount, 0);
long fileSampleCount = Header.LengthInSamples;
long fileSampleOffset = -delaySampleCount;
int chunkSampleCount = nfft * bitmapWidth;
int chunkCount = (int)Math.Ceiling((double)(fileSampleCount + delaySampleCount) / chunkSampleCount);
double[] chunkSamples = new double[chunkSampleCount];
_data = new byte[chunkSampleCount * Header.BlockAlign];
_stream.Seek(Header.DataStartPosition, SeekOrigin.Begin);
// for negative delays, skip samples at the beginning
if (fileSampleOffset > 0)
{
samples.Add(0);
if (samples.Count == sampleSize)
{
var samplesAsReal = new double[sampleSize];
for (int k = 0; k < sampleSize; k++)
samplesAsReal[k] = 0;
var bmp = DrawSpectrogram(nfft, samplesAsReal, f, palette);
bmp.Save(Path.Combine(spectrogramDirectory, count + ".gif"), System.Drawing.Imaging.ImageFormat.Gif);
bitmaps.Add(bmp);
samples = new List<int>();
count++;
}
_stream.Seek(fileSampleOffset * Header.BlockAlign, SeekOrigin.Current);
}
// load data in smaller parts
_data = new byte[Header.BytesPerSecond];
_stream.Position = Header.DataStartPosition;
int bytesRead = _stream.Read(_data, 0, _data.Length);
while (bytesRead == Header.BytesPerSecond)
for (int iChunk = 0; iChunk < chunkCount; iChunk++)
{
while (index < Header.BytesPerSecond)
// calculate padding at the beginning (for positive delays)
int startPaddingSampleCount = 0;
if (fileSampleOffset < 0)
{
int value = 0;
for (int channelNumber = 0; channelNumber < Header.NumberOfChannels; channelNumber++)
{
value += readSampleDataValue.Invoke(ref index);
}
value = value / Header.NumberOfChannels;
if (value < DataMinValue)
DataMinValue = value;
if (value > DataMaxValue)
DataMaxValue = value;
samples.Add(value);
totalSamples++;
startPaddingSampleCount = (int)Math.Min(-fileSampleOffset, chunkSampleCount);
fileSampleOffset += startPaddingSampleCount;
}
if (samples.Count == sampleSize)
// calculate how many samples to read from the file
long fileSamplesRemaining = fileSampleCount - Math.Max(fileSampleOffset, 0);
int fileReadSampleCount = (int)Math.Min(fileSamplesRemaining, chunkSampleCount - startPaddingSampleCount);
// calculate padding at the end (when the data isn't an even multiple of our chunk size)
int endPaddingSampleCount = chunkSampleCount - startPaddingSampleCount - fileReadSampleCount;
int chunkSampleOffset = 0;
// add padding at the beginning
if (startPaddingSampleCount > 0)
{
Array.Clear(chunkSamples, chunkSampleOffset, startPaddingSampleCount);
chunkSampleOffset += startPaddingSampleCount;
}
// read samples from the file
if (fileReadSampleCount > 0)
{
int fileReadByteCount = fileReadSampleCount * Header.BlockAlign;
_stream.Read(_data, 0, fileReadByteCount);
fileSampleOffset += fileReadSampleCount;
int dataByteOffset = 0;
while (dataByteOffset < fileReadByteCount)
{
var samplesAsReal = new double[sampleSize];
for (int k = 0; k < sampleSize; k++)
samplesAsReal[k] = samples[k] / divider;
var bmp = DrawSpectrogram(nfft, samplesAsReal, f, palette);
bmp.Save(Path.Combine(spectrogramDirectory, count + ".gif"), System.Drawing.Imaging.ImageFormat.Gif);
bitmaps.Add(bmp);
samples = new List<int>();
count++;
int value = 0;
for (int iChannel = 0; iChannel < Header.NumberOfChannels; iChannel++)
{
value += readSampleDataValue(ref dataByteOffset);
}
chunkSamples[chunkSampleOffset] = value * sampleScale;
chunkSampleOffset += 1;
}
}
bytesRead = _stream.Read(_data, 0, _data.Length);
index = 0;
}
if (samples.Count > 0)
{
var samplesAsReal = new double[sampleSize];
for (int k = 0; k < sampleSize && k < samples.Count; k++)
samplesAsReal[k] = samples[k] / divider;
var bmp = DrawSpectrogram(nfft, samplesAsReal, f, palette);
bmp.Save(Path.Combine(spectrogramDirectory, count + ".gif"), System.Drawing.Imaging.ImageFormat.Gif);
// add padding at the end
if (endPaddingSampleCount > 0)
{
Array.Clear(chunkSamples, chunkSampleOffset, endPaddingSampleCount);
chunkSampleOffset += endPaddingSampleCount;
}
// generate spectrogram for this chunk
Bitmap bmp = drawer.Draw(chunkSamples);
bmp.Save(Path.Combine(spectrogramDirectory, iChunk + ".gif"), System.Drawing.Imaging.ImageFormat.Gif);
bitmaps.Add(bmp);
}
var doc = new XmlDocument();
doc.LoadXml("<SpectrogramInfo><SampleDuration/><TotalDuration/><AudioFormat /><AudioFormat /><ChunkId /><SecondsPerImage /><ImageWidth /><NFFT /></SpectrogramInfo>");
double sampleDuration = Header.LengthInSeconds / (totalSamples / Convert.ToDouble(nfft));
doc.DocumentElement.SelectSingleNode("SampleDuration").InnerText = sampleDuration.ToString(CultureInfo.InvariantCulture);
doc.DocumentElement.SelectSingleNode("SampleDuration").InnerText = ((double)nfft / Header.SampleRate).ToString(CultureInfo.InvariantCulture);
doc.DocumentElement.SelectSingleNode("TotalDuration").InnerText = Header.LengthInSeconds.ToString(CultureInfo.InvariantCulture);
doc.DocumentElement.SelectSingleNode("AudioFormat").InnerText = Header.AudioFormat.ToString(CultureInfo.InvariantCulture);
doc.DocumentElement.SelectSingleNode("ChunkId").InnerText = Header.ChunkId.ToString(CultureInfo.InvariantCulture);
doc.DocumentElement.SelectSingleNode("SecondsPerImage").InnerText = ((double)(sampleSize / (double)Header.SampleRate)).ToString(CultureInfo.InvariantCulture);
doc.DocumentElement.SelectSingleNode("SecondsPerImage").InnerText = ((double)chunkSampleCount / Header.SampleRate).ToString(CultureInfo.InvariantCulture);
doc.DocumentElement.SelectSingleNode("ImageWidth").InnerText = bitmapWidth.ToString(CultureInfo.InvariantCulture);
doc.DocumentElement.SelectSingleNode("NFFT").InnerText = nfft.ToString(CultureInfo.InvariantCulture);
doc.Save(Path.Combine(spectrogramDirectory, "Info.xml"));
@ -508,131 +501,202 @@ namespace Nikse.SubtitleEdit.Core
return bitmaps;
}
private static Bitmap DrawSpectrogram(int nfft, double[] samples, Fourier f, Color[] palette)
private class SpectrogramDrawer
{
const int overlap = 0;
int numSamples = samples.Length;
int colIncrement = nfft * (1 - overlap);
private const double raisedCosineWindowScale = 0.5;
int numcols = numSamples / colIncrement;
// make sure we don't step beyond the end of the recording
while ((numcols - 1) * colIncrement + nfft > numSamples)
numcols--;
private int _nfft;
private MagnitudeToIndexMapper _mapper;
private RealFFT _fft;
private FastBitmap.PixelData[] _palette;
private double[] _segment;
private double[] _window;
private double[] _magnitude;
double[] real = new double[nfft];
double[] imag = new double[nfft];
double[] magnitude = new double[nfft / 2];
var bmp = new Bitmap(numcols, nfft / 2);
for (int col = 0; col <= numcols - 1; col++)
public SpectrogramDrawer(int nfft)
{
// read a segment of the recorded signal
for (int c = 0; c <= nfft - 1; c++)
_nfft = nfft;
_mapper = new MagnitudeToIndexMapper(100.0, 255);
_fft = new RealFFT(nfft);
_palette = GeneratePalette(nfft);
_segment = new double[nfft];
_window = CreateRaisedCosineWindow(nfft);
_magnitude = new double[nfft / 2];
double scaleCorrection = 1.0 / (raisedCosineWindowScale * _fft.ForwardScaleFactor);
for (int i = 0; i < _window.Length; i++)
{
imag[c] = 0;
real[c] = samples[col * colIncrement + c] * Fourier.Hanning(nfft, c);
}
// transform to the frequency domain
f.FourierTransform(real, imag);
// and compute the magnitude spectrum
f.MagnitudeSpectrum(real, imag, Fourier.W0Hanning, magnitude);
// Draw
for (int newY = 0; newY < nfft / 2 - 1; newY++)
{
int colorIndex = MapToPixelIndex(magnitude[newY], 100, 255);
bmp.SetPixel(col, (nfft / 2 - 1) - newY, palette[colorIndex]);
_window[i] *= scaleCorrection;
}
}
public Bitmap Draw(double[] samples)
{
int width = samples.Length / _nfft;
int height = _nfft / 2;
var bmp = new FastBitmap(new Bitmap(width, height));
bmp.LockImage();
for (int x = 0; x < width; x++)
{
// read a segment of the recorded signal
for (int i = 0; i < _nfft; i++)
{
_segment[i] = samples[x * _nfft + i] * _window[i];
}
// transform to the frequency domain
_fft.ComputeForward(_segment);
// compute the magnitude of the spectrum
MagnitudeSpectrum(_segment, _magnitude);
// draw
for (int y = 0; y < height; y++)
{
int colorIndex = _mapper.Map(_magnitude[y]);
bmp.SetPixel(x, height - y - 1, _palette[colorIndex]);
}
}
bmp.UnlockImage();
return bmp.GetBitmap();
}
private static double[] CreateRaisedCosineWindow(int n)
{
double twoPiOverN = Math.PI * 2.0 / n;
double[] dst = new double[n];
for (int i = 0; i < n; i++)
dst[i] = 0.5 * (1.0 - Math.Cos(twoPiOverN * i));
return dst;
}
private static void MagnitudeSpectrum(double[] segment, double[] magnitude)
{
magnitude[0] = Math.Sqrt(SquareSum(segment[0], segment[1]));
for (int i = 2; i < segment.Length; i += 2)
magnitude[i / 2] = Math.Sqrt(SquareSum(segment[i], segment[i + 1]) * 2.0);
}
private static double SquareSum(double a, double b)
{
return a * a + b * b;
}
private static FastBitmap.PixelData[] GeneratePalette(int nfft)
{
var palette = new FastBitmap.PixelData[nfft];
if (Configuration.Settings.VideoControls.SpectrogramAppearance == "Classic")
{
for (int colorIndex = 0; colorIndex < nfft; colorIndex++)
palette[colorIndex] = new FastBitmap.PixelData(PaletteValue(colorIndex, nfft));
}
else
{
var list = SmoothColors(0, 0, 0, Configuration.Settings.VideoControls.WaveformColor.R,
Configuration.Settings.VideoControls.WaveformColor.G,
Configuration.Settings.VideoControls.WaveformColor.B, nfft);
for (int i = 0; i < nfft; i++)
palette[i] = new FastBitmap.PixelData(list[i]);
}
return palette;
}
private static Color PaletteValue(int x, int range)
{
double g;
double r;
double b;
double r4 = range / 4.0;
const double u = 255;
if (x < r4)
{
b = x / r4;
g = 0;
r = 0;
}
else if (x < 2 * r4)
{
b = (1 - (x - r4) / r4);
g = 1 - b;
r = 0;
}
else if (x < 3 * r4)
{
b = 0;
g = (2 - (x - r4) / r4);
r = 1 - g;
}
else
{
b = (x - 3 * r4) / r4;
g = 0;
r = 1 - b;
}
r = ((int)(Math.Sqrt(r) * u)) & 0xff;
g = ((int)(Math.Sqrt(g) * u)) & 0xff;
b = ((int)(Math.Sqrt(b) * u)) & 0xff;
return Color.FromArgb((int)r, (int)g, (int)b);
}
private static List<Color> SmoothColors(int fromR, int fromG, int fromB, int toR, int toG, int toB, int count)
{
while (toR < 255 && toG < 255 && toB < 255)
{
toR++;
toG++;
toB++;
}
var list = new List<Color>();
double r = fromR;
double g = fromG;
double b = fromB;
double diffR = (toR - fromR) / (double)count;
double diffG = (toG - fromG) / (double)count;
double diffB = (toB - fromB) / (double)count;
for (int i = 0; i < count; i++)
{
list.Add(Color.FromArgb((int)r, (int)g, (int)b));
r += diffR;
g += diffG;
b += diffB;
}
return list;
}
/// Maps magnitudes in the range [-decibelRange .. 0] dB to palette index values in the range [0 .. indexRange-1]
private class MagnitudeToIndexMapper
{
private readonly double _minMagnitude;
private readonly double _multiplier;
private readonly double _addend;
public MagnitudeToIndexMapper(double decibelRange, int indexRange)
{
double mappingScale = indexRange / decibelRange;
_minMagnitude = Math.Pow(10.0, -decibelRange / 20.0);
_multiplier = 20.0 * mappingScale;
_addend = decibelRange * mappingScale;
}
public int Map(double magnitude)
{
return magnitude >= _minMagnitude ? (int)(Math.Log10(magnitude) * _multiplier + _addend) : 0;
}
// Less optimized but readable version of the above
public static int Map(double magnitude, double decibelRange, int indexRange)
{
if (magnitude == 0) return 0;
double decibelLevel = 20.0 * Math.Log10(magnitude);
return decibelLevel >= -decibelRange ? (int)(indexRange * (decibelLevel + decibelRange) / decibelRange) : 0;
}
}
return bmp;
}
public static Color PaletteValue(int x, int range)
{
double g;
double r;
double b;
double r4 = range / 4.0;
const double u = 255;
if (x < r4)
{
b = x / r4;
g = 0;
r = 0;
}
else if (x < 2 * r4)
{
b = (1 - (x - r4) / r4);
g = 1 - b;
r = 0;
}
else if (x < 3 * r4)
{
b = 0;
g = (2 - (x - r4) / r4);
r = 1 - g;
}
else
{
b = (x - 3 * r4) / r4;
g = 0;
r = 1 - b;
}
r = ((int)(Math.Sqrt(r) * u)) & 0xff;
g = ((int)(Math.Sqrt(g) * u)) & 0xff;
b = ((int)(Math.Sqrt(b) * u)) & 0xff;
return Color.FromArgb((int)r, (int)g, (int)b);
}
/// <summary>
/// Maps magnitudes in the range [-rangedB .. 0] dB to palette index values in the range [0 .. rangeIndex-1]
/// and computes and returns the index value which corresponds to passed-in magnitude
/// </summary>
private static int MapToPixelIndex(double magnitude, double rangedB, int rangeIndex)
{
const double log10 = 2.30258509299405;
if (magnitude == 0)
return 0;
double levelIndB = 20 * Math.Log(magnitude) / log10;
if (levelIndB < -rangedB)
return 0;
return (int)(rangeIndex * (levelIndB + rangedB) / rangedB);
}
private static List<Color> SmoothColors(int fromR, int fromG, int fromB, int toR, int toG, int toB, int count)
{
while (toR < 255 && toG < 255 && toB < 255)
{
toR++;
toG++;
toB++;
}
var list = new List<Color>();
double r = fromR;
double g = fromG;
double b = fromB;
double diffR = (toR - fromR) / (double)count;
double diffG = (toG - fromG) / (double)count;
double diffB = (toB - fromB) / (double)count;
for (int i = 0; i < count; i++)
{
list.Add(Color.FromArgb((int)r, (int)g, (int)b));
r += diffR;
g += diffG;
b += diffB;
}
return list;
}
}
}