1 /* 2 * Copyright (c) 2013 The WebRTC project authors. All Rights Reserved. 3 * 4 * Use of this source code is governed by a BSD-style license 5 * that can be found in the LICENSE file in the root of the source 6 * tree. An additional intellectual property rights grant can be found 7 * in the file PATENTS. All contributing project authors may 8 * be found in the AUTHORS file in the root of the source tree. 9 */ 10 11 #include "dl/sp/src/test/gensig.h" 12 13 #include <stdio.h> 14 #include <math.h> 15 #include <stdlib.h> 16 17 #define MAX_REAL_SIGNAL_TYPE 3 18 #define MAX_SIGNAL_TYPE 4 19 20 int MaxSignalType(int real_only) { 21 return real_only ? MAX_REAL_SIGNAL_TYPE : MAX_SIGNAL_TYPE; 22 } 23 24 /* 25 * Generate a test signal and compute the theoretical FFT. 26 * 27 * The test signal is specified by |signal_type|, and the test signal 28 * is saved in |x| with the corresponding FFT in |fft|. The size of 29 * the test signal is |size|. |signalValue| is desired the amplitude 30 * of the test signal. 31 * 32 * If |real_only| is true, then the test signal is assumed to be real 33 * instead of complex, which is the default. This is only applicable 34 * for a |signal_type| of 0 or 3; the other signals are already real-valued. 35 */ 36 void GenerateTestSignalAndFFT(struct ComplexFloat* x, 37 struct ComplexFloat* fft, 38 int size, 39 int signal_type, 40 float signal_value, 41 int real_only) { 42 int k; 43 44 switch (signal_type) { 45 case 0: 46 /* 47 * Signal is a constant signal_value + i*signal_value (or just 48 * signal_value if real_only is true.) 49 */ 50 for (k = 0; k < size; ++k) { 51 x[k].Re = signal_value; 52 x[k].Im = real_only ? 0 : signal_value; 53 } 54 55 fft[0].Re = signal_value * size; 56 fft[0].Im = real_only ? 0 : signal_value * size; 57 58 for (k = 1; k < size; ++k) { 59 fft[k].Re = fft[k].Im = 0; 60 } 61 break; 62 case 1: 63 /* 64 * A real-valued ramp 65 */ 66 { 67 double factor = signal_value / (float) size; 68 double omega = 2 * M_PI / size; 69 70 for (k = 0; k < size; ++k) { 71 x[k].Re = ((k + 1)*factor); 72 x[k].Im = 0; 73 } 74 75 fft[0].Re = factor * size * (size + 1) / 2; 76 fft[0].Im = 0; 77 for (k = 1; k < size; ++k) { 78 double phase; 79 phase = omega * k; 80 fft[k].Re = factor * -size / 2; 81 fft[k].Im = factor * size / 2 * (sin(phase) / (1 - cos(phase))); 82 } 83 84 /* 85 * Remove any roundoff for k = N/2 since sin(2*pi/N*N/2) = 0. 86 */ 87 fft[size / 2].Im = 0; 88 } 89 break; 90 case 2: 91 /* 92 * Pure real-valued sine wave, one cycle. 93 */ 94 { 95 double omega = 2 * M_PI / size; 96 97 for (k = 0; k < size; ++k) { 98 x[k].Re = signal_value * sin(omega * k); 99 x[k].Im = 0; 100 } 101 102 /* 103 * Remove any roundoff for k = N/2 since sin(2*pi/N*N/2) = 0. 104 */ 105 x[size / 2 ].Re = 0; 106 107 for (k = 0; k < size; ++k) { 108 fft[k].Re = 0; 109 fft[k].Im = 0; 110 } 111 112 /* 113 * When size == 2, x[k] is identically zero, so the FFT is also zero. 114 */ 115 if (size != 2) { 116 fft[1].Im = -signal_value * (size / 2); 117 fft[size - 1].Im = signal_value * (size / 2); 118 } 119 } 120 break; 121 case 3: 122 /* 123 * The source signal is x[k] = 0 except x[1] = x[size-1] = 124 * -i*signal_value. The transform is one period of a pure real 125 * (negative) sine wave. Only defined when real_only is false. 126 */ 127 if (!real_only) { 128 double omega = 2 * M_PI / size; 129 for (k = 0; k < size; ++k) { 130 x[k].Re = 0; 131 x[k].Im = 0; 132 } 133 x[1].Im = -signal_value; 134 x[size-1].Im = signal_value; 135 136 if (size == 2) { 137 fft[0].Re = 0; 138 fft[0].Im = signal_value; 139 fft[1].Re = 0; 140 fft[1].Im = -signal_value; 141 } else { 142 for (k = 0; k < size; ++k) { 143 fft[k].Re = -2 * signal_value * sin(omega * k); 144 fft[k].Im = 0; 145 } 146 147 /* 148 * Remove any roundoff for k = N/2 since sin(2*pi/N*N/2) = 0. 149 */ 150 fft[size / 2].Re = 0; 151 } 152 break; 153 } 154 /* Fall through if real_only */ 155 case MAX_SIGNAL_TYPE: 156 default: 157 fprintf(stderr, "invalid signal type: %d\n", signal_type); 158 exit(1); 159 } 160 } 161