想当然了,用递归实现DCT,没想到DCT有4个分支需要递归下去,这样的规模非但无法快速实现,反而由于本身时间复杂度没有多大减少加上递归开销等等比慢速实现往往还慢。
这个代码片段将由于清洁需要从QSharp中删除而保留在这里,对其分析将在代码之后有空时进行。过两天想想是不是能用动态规划或备忘录来改进这个算法。
/// <summary>
/// Type-IV DCT implemented using recursive method
/// Warning: This method is mathematically crippled even slower than the normal one and thereby not suitable for fast caluculation
/// </summary>
/// <param name="input">The input sequence</param>
/// <param name="output">The output sequence</param>
public static void FctIvRecursive(this IList<double> input, IList<double> output)
{
var n = new DftMethods.FactorAidedInteger(input.Count);
input.FctIvRecursiveCore(n, 0, 1, output);
}
private static void DctIvFCore(this IList<double> input, int start, int skip, double a, IList<double> output)
{
var ni = input.Count;
var n = output.Count;
for (var k = 0; k < n; k++)
{
var d = 0.0;
var j = 0;
for (var i = start; i < ni; i += skip, j++)
{
var v = input[i] * Math.Cos(a * (k + 0.5) * (j + 0.5));
if (j % 2 != 0)
{
v = -v;
}
d += v;
}
output[k] = d;
}
}
private static void DstIvCore(this IList<double> input, int start, int skip, double a, IList<double> output)
{
var ni = input.Count;
var n = output.Count;
for (var k = 0; k < n; k++)
{
var d = 0.0;
var j = 0;
for (var i = start; i < ni; i += skip, j++)
{
d += input[i]*Math.Sin(a*(k + 0.5)*(j + 0.5));
}
output[k] = d;
}
}
private static void DstIvFCore(this IList<double> input, int start, int skip, double a, IList<double> output)
{
var ni = input.Count;
var n = output.Count;
for (var k = 0; k < n; k++)
{
var d = 0.0;
var j = 0;
for (var i = start; i < ni; i += skip, j++)
{
var v = input[i] * Math.Sin(a * (k + 0.5) * (j + 0.5));
if (j % 2 != 0)
{
v = -v;
}
d += v;
}
output[k] = d;
}
}
private static void FctIvRecursiveCore(this IList<double> input,
DftMethods.FactorAidedInteger n, int start, int skip,
IList<double> output)
{
var n0 = n.N;
var n1 = n.Decompose();
var n2 = n.N;
var a = Math.PI/n0;
if (n1 > 1 && n2 > 1)
{
for (var k = 0; k < n0; k++)
{
output[k] = 0;
}
var b = 0.25*a;
for (var r = 0; r < n1; r++)
{
int newStart, newSkip;
DftMethods.GetSubsequenceIndexGenerator(start, skip, r, n1, out newStart, out newSkip);
// N2-point DCT and DST
var tmpDct = new double[n2];
var tmpDst = new double[n2];
var tmpDctF = new double[n2];
var tmpDstF = new double[n2];
input.FctIvRecursiveCore(n, newStart, newSkip, tmpDct);
input.FstIvRecursiveCore(n, newStart, newSkip, tmpDst);
input.FctIvFRecursiveCore(n, newStart, newSkip, tmpDctF);
input.FstIvFRecursiveCore(n, newStart, newSkip, tmpDstF);
var c = b*(2*r + 1 - n1);
for (var k = 0; k < n0; k++)
{
var k2 = k%n2;
var l = k/n2;
var phi = c*(2*k + 1);
var lMod4 = l%4;
switch (lMod4)
{
case 0:
output[k] += tmpDct[k2] * Math.Cos(phi) - tmpDst[k2] * Math.Sin(phi);
break;
case 2:
output[k] += -tmpDct[k2] * Math.Cos(phi) + tmpDst[k2] * Math.Sin(phi);
break;
case 1:
output[k] += -tmpDstF[k2] * Math.Cos(phi) - tmpDctF[k2] * Math.Sin(phi);
break;
case 3:
output[k] += tmpDstF[k2] * Math.Cos(phi) + tmpDctF[k2] * Math.Sin(phi);
break;
}
}
}
}
else
{
input.DctIvCore(start, skip, a, output);
}
n.Restore(n1);
}
private static void FstIvRecursiveCore(this IList<double> input,
DftMethods.FactorAidedInteger n, int start, int skip,
IList<double> output)
{
var n0 = n.N;
var n1 = n.Decompose();
var n2 = n.N;
var a = Math.PI/n0;
if (n1 > 1 && n2 > 1)
{
for (var k = 0; k < n0; k++)
{
output[k] = 0;
}
var b = 0.25*a;
for (var r = 0; r < n1; r++)
{
int newStart, newSkip;
DftMethods.GetSubsequenceIndexGenerator(start, skip, r, n1, out newStart, out newSkip);
// N2-point DCT and DST
var tmpDct = new double[n2];
var tmpDst = new double[n2];
var tmpDctF = new double[n2];
var tmpDstF = new double[n2];
input.FctIvRecursiveCore(n, newStart, newSkip, tmpDct);
input.FstIvRecursiveCore(n, newStart, newSkip, tmpDst);
input.FctIvFRecursiveCore(n, newStart, newSkip, tmpDctF);
input.FstIvFRecursiveCore(n, newStart, newSkip, tmpDstF);
var c = b*(2*r + 1 - n1);
for (var k = 0; k < n0; k++)
{
var k2 = k%n2;
var l = k/n2;
var phi = c*(2*k + 1);
var lMod4 = l % 4;
switch (lMod4)
{
case 0:
output[k] += tmpDst[k2] * Math.Cos(phi) + tmpDct[k2] * Math.Sin(phi);
break;
case 2:
output[k] += -tmpDst[k2] * Math.Cos(phi) - tmpDct[k2] * Math.Sin(phi);
break;
case 1:
output[k] += tmpDctF[k2] * Math.Cos(phi) - tmpDstF[k2] * Math.Sin(phi);
break;
case 3:
output[k] += -tmpDctF[k2] * Math.Cos(phi) + tmpDstF[k2] * Math.Sin(phi);
break;
}
}
}
}
else
{
input.DstIvCore(start, skip, a, output);
}
n.Restore(n1);
}
private static void FctIvFRecursiveCore(this IList<double> input,
DftMethods.FactorAidedInteger n, int start, int skip,
IList<double> output)
{
var n0 = n.N;
var n1 = n.Decompose();
var n2 = n.N;
var a = Math.PI / n0;
if (n1 > 1 && n2 > 1)
{
for (var k = 0; k < n0; k++)
{
output[k] = 0;
}
var b = 0.25 * a;
for (var r = 0; r < n1; r++)
{
int newStart, newSkip;
DftMethods.GetSubsequenceIndexGenerator(start, skip, r, n1, out newStart, out newSkip);
// N2-point DCT and DST
var tmpDct1 = new double[n2];
var tmpDst1 = new double[n2];
var tmpDct2 = new double[n2];
var tmpDst2 = new double[n2];
if (n1%2 == 0)
{
input.FctIvRecursiveCore(n, newStart, newSkip, tmpDct1);
input.FstIvRecursiveCore(n, newStart, newSkip, tmpDst1);
input.FctIvFRecursiveCore(n, newStart, newSkip, tmpDct2);
input.FstIvFRecursiveCore(n, newStart, newSkip, tmpDst2);
}
else
{
input.FctIvRecursiveCore(n, newStart, newSkip, tmpDct2);
input.FstIvRecursiveCore(n, newStart, newSkip, tmpDst2);
input.FctIvFRecursiveCore(n, newStart, newSkip, tmpDct1);
input.FstIvFRecursiveCore(n, newStart, newSkip, tmpDst1);
}
var c = b*(2*r + 1 - n1);
var rIsOdd = r%2 != 0;
for (var k = 0; k < n0; k++)
{
var k2 = k % n2;
var l = k / n2;
var phi = c * (2 * k + 1);
var lMod4 = l % 4;
var d = 0.0;
switch (lMod4)
{
case 0:
d = tmpDct1[k2] * Math.Cos(phi) - tmpDst1[k2] * Math.Sin(phi);
break;
case 2:
d = -tmpDct1[k2] * Math.Cos(phi) + tmpDst1[k2] * Math.Sin(phi);
break;
case 1:
d = -tmpDst2[k2] * Math.Cos(phi) - tmpDct2[k2] * Math.Sin(phi);
break;
case 3:
d = tmpDst2[k2] * Math.Cos(phi) + tmpDct2[k2] * Math.Sin(phi);
break;
}
if (rIsOdd)
{
d = -d;
}
output[k] += d;
}
}
}
else
{
input.DctIvFCore(start, skip, a, output);
}
n.Restore(n1);
}
private static void FstIvFRecursiveCore(this IList<double> input,
DftMethods.FactorAidedInteger n, int start, int skip,
IList<double> output)
{
var n0 = n.N;
var n1 = n.Decompose();
var n2 = n.N;
var a = Math.PI / n0;
if (n1 > 1 && n2 > 1)
{
for (var k = 0; k < n0; k++)
{
output[k] = 0;
}
var b = 0.25 * a;
for (var r = 0; r < n1; r++)
{
int newStart, newSkip;
DftMethods.GetSubsequenceIndexGenerator(start, skip, r, n1, out newStart, out newSkip);
// N2-point DCT and DST
var tmpDct1 = new double[n2];
var tmpDst1 = new double[n2];
var tmpDct2 = new double[n2];
var tmpDst2 = new double[n2];
if (n1%2 == 0)
{
input.FctIvRecursiveCore(n, newStart, newSkip, tmpDct1);
input.FstIvRecursiveCore(n, newStart, newSkip, tmpDst1);
input.FctIvFRecursiveCore(n, newStart, newSkip, tmpDct2);
input.FstIvFRecursiveCore(n, newStart, newSkip, tmpDst2);
}
else
{
input.FctIvRecursiveCore(n, newStart, newSkip, tmpDct2);
input.FstIvRecursiveCore(n, newStart, newSkip, tmpDst2);
input.FctIvFRecursiveCore(n, newStart, newSkip, tmpDct1);
input.FstIvFRecursiveCore(n, newStart, newSkip, tmpDst1);
}
var c = b * (2 * r + 1 - n1);
var rIsOdd = r%2 != 0;
for (var k = 0; k < n0; k++)
{
var k2 = k % n2;
var l = k / n2;
var phi = c * (2 * k + 1);
var lMod4 = l % 4;
var d = 0.0;
switch (lMod4)
{
case 0:
d = tmpDst1[k2] * Math.Cos(phi) + tmpDct1[k2] * Math.Sin(phi);
break;
case 2:
d = -tmpDst1[k2] * Math.Cos(phi) - tmpDct1[k2] * Math.Sin(phi);
break;
case 1:
d = tmpDct2[k2] * Math.Cos(phi) - tmpDst2[k2] * Math.Sin(phi);
break;
case 3:
d = tmpDct2[k2] * Math.Cos(phi) - tmpDst2[k2] * Math.Sin(phi);
break;
}
if (rIsOdd)
{
d = -d;
}
output[k] += d;
}
}
}
else
{
input.DstIvFCore(start, skip, a, output);
}
n.Restore(n1);
}