1. 時域與頻域
時域（時間域）——自變量是時間,即橫軸是時間,縱軸是信號的變化。其動態信號x（t）是描述信號在不同時刻取值的函數；
頻域（頻率域）——自變量是頻率,即橫軸是頻率,縱軸是該頻率信號的幅度,也就是通常說的頻譜圖。
正弦信號： x(t) = A*sin(w*t+p); 其中，A為幅值；頻率為f;W為角頻率=2pi*f；P為初始相位； 周期T=1/f；
注：實際物理頻率f:表示AD采集物理信號的頻率；fs為采樣頻率，由奈奎斯特采樣定理可以知道，fs必須≥信號最高頻率的2倍才不會發生信號混疊，因此fs能采樣到的信號最高頻率為fs/2。
角頻率是物理頻率的2*pi倍，這個也稱模擬頻率；（由於一個信號周期（如交流電）是360度，即2pi。故角頻率就是轉了多少個2pi。設置角頻率純粹為了便於計算）。
歸一化頻率是將實際物理頻率按fs歸一化之后的結果，最高的信號頻率為fs/2對應歸一化頻率0.5，這也就是為什么在matlab的fdtool工具中歸一化頻率為什么最大只到0.5的原因。
圓周頻率是歸一化頻率的2*pi倍，這個也稱數字頻率。也就是歸一化的角頻率。
如：x(t) = 5*sin(2*pi*6*t+0/180*pi);
幅值A為5; 頻率f為6Hz； 角頻率W=2*pi*6(rad/s);初始相位為0；周期T為1/6；
圖1是正弦波的時域圖，示出了振幅與時間的關系。
- 在時域圖中，橫軸是時間，縱軸是振幅。
- 時域圖顯示振幅隨時間的變化，可以看出峰值振幅為5V，可以算出頻率f=6 Hz。
圖2是圖1中正弦波的頻域圖
- 在頻域圖中，橫軸是頻率，縱軸是峰值振幅。
- 頻域圖僅僅示出峰值振幅與頻率，而不顯示振幅隨時間的變化。
- 從頻域圖可以看出，正弦波的頻率為6Hz，這個6Hz的正弦波的峰值振幅為5V
- 頻域圖的優點是，從頻域圖中，可以一眼看出正弦波的頻率和峰值振幅
- 整個正弦波在頻域圖上只是一個立柱
- 立柱的位置顯示了正弦波的頻率
- 立柱的高度顯示了正弦波的峰值振幅
2. 傅里葉變換
傅里葉變換是把時域信號轉到頻域，進行頻譜分析、濾波處理的關鍵一步，掌握傅里葉變換的過程、算法原理及時域與頻域之間的聯系，有助於后續進行數字濾波器的設計。
圖3 時域圖像與頻域圖像
傅里葉變換的本質、算法原理參見我的另一篇博客：https://www.cnblogs.com/zhaopengpeng/p/14476148.html
3. 頻域分析
頻域分析使用的是傅里葉變換。根據傅里葉定理，任何連續測量的時序或信號f(t)，都可以表示為不同頻率的正弦波信號的無限疊加.即
混合信號的頻譜分析，詳見可見我的另一博客：https://www.cnblogs.com/zhaopengpeng/p/15224552.html
4.信號的傅里葉變換FFT與實際的采樣頻率間的關系
- 當已知數據實際的采樣頻率（比如Fs=1200Hz）時，可確定數據對應的實際頻率范圍為0-1200Hz，根據奈奎斯特定理可知，當頻率f=Fs/2范圍[0,600]內的信號才是被采樣到的有效信號；
- 做N個點的FFT，表示在時域上對原來的信號取了N個點來做頻譜分析，N點FFT變換后的結果仍為N個點（這N個點是復數形式，極坐標系下，與實部的夾角相當於頻率，向量的長度相當於幅值）；
- FFT變換后的N個點，每個點對應的實際頻率是多少，必須知道原來信號的采樣頻率Fs，才能計算第k個點的實際頻率，即f(k)=k*(Fs/N), k=0,1,2,...,N; 其中的Fs/N就是頻率分辨率（=采樣頻率/FFT處理的數據點個數）；FFT時，也可以不用原來信號的采樣頻率，重新采用新的采樣頻率對數據進行采樣，但是這就無法保證FFT變換后的N個點對應的實際頻率，與重采樣后計算的頻率保持一致了；
- 當Fs=1200Hz，N=256個點做FFT之后，因為頻率分辨率為4.6875Hz，如果原來的信號在50Hz或100Hz有分量時，在頻譜上是看不見的，就必須取越多的點數來做FFT，N就越大，在時域上就必須取更長的信號樣本做分析。
5. 信號的頻域分析及合成代碼（Matlab與C++實現）
Matlab代碼：
1 clear all; close all; 2 3 % 外界振動信號的模擬 4 5 A0 = 745200; 6 A1 = 5636; %頻率F1信號的幅度 7 A2 = 4950; %頻率F2信號的幅度 8 A3 = 2422; %頻率F2信號的幅度 9 A4 = 2372; %頻率F2信號的幅度 10 A5 = 2268; %頻率F2信號的幅度 11 A6 = 1466; %頻率F1信號的幅度 12 A7 = 1394; %頻率F2信號的幅度 13 A8 = 1177; %頻率F2信號的幅度 14 A9 = 1108; %頻率F2信號的幅度 15 A10 = 938.2; %頻率F2信號的幅度 16 A11 = 929.1; %頻率F2信號的幅度 17 18 F1 = 98.44; %信號1頻率(Hz) 19 F2 = 51.56; %信號2頻率(Hz) 20 F3 = 103.1; 21 F4 = 150; 22 F5 = 46.88; 23 F6 = 93.75; %信號1頻率(Hz) 24 F7 = 56.25; %信號2頻率(Hz) 25 F8 = 201.6; 26 F9 = 284.4; 27 F10 = 107.8; 28 F11 = 398.4; 29 30 P1 = 31.87; %信號1相位(度) 31 P2 = 84.84; %信號2相位(度) 32 P3 = -151.8; 33 P4 = 17.4; 34 P5 = -90.85; 35 P6 = 35.24; %信號1相位(度) 36 P7 = 75.83 %信號2相位(度) 37 P8 = 135.2; 38 P9 = 166.7; 39 P10 = -157.6; 40 P11 = 79.09; 41 42 Fs = 1200; %采樣頻率(Hz) 43 44 filename = 'E:\BYDA\trunk\build\bin\scale_data_static\202108141200\scale_data_static_cov\sinc4_1200sps_0kg.dat'; 45 [data] = textread(filename, '%n'); 46 % 47 N =length(data); 48 % N = 256; %采樣點數，請注意時域采樣N點，FFT時也是N點！ 49 t = [0 : 1/Fs: N/Fs]; %采樣時刻，注意不是序列是真正的時間！t = [0:N]*Ts 50 51 % 生成信號 52 signal1=A1*cos(2*pi*F1*t+pi*P1/180); % 周期T1 = 2*pi/(2*pi*F1) = 0.01s，采樣周期Ts = 0.00019531 < (T1)/2 53 signal2=A2*cos(2*pi*F2*t+pi*P2/180); % 周期T2 = 2*pi/(2*pi*F2) = 0.001s，采樣周期Ts = 0.00019531 < (T2)/2 54 signal3=A3*cos(2*pi*F3*t+pi*P3/180); 55 signal4=A4*cos(2*pi*F4*t+pi*P4/180); 56 signal5=A5*cos(2*pi*F5*t+pi*P5/180); 57 58 signal6=A6*cos(2*pi*F6*t+pi*P6/180); % 周期T1 = 2*pi/(2*pi*F1) = 0.01s，采樣周期Ts = 0.00019531 < (T1)/2 59 signal7=A7*cos(2*pi*F7*t+pi*P7/180); % 周期T2 = 2*pi/(2*pi*F2) = 0.001s，采樣周期Ts = 0.00019531 < (T2)/2 60 signal8=A8*cos(2*pi*F8*t+pi*P8/180); 61 signal9=A9*cos(2*pi*F9*t+pi*P9/180); 62 signal10=A10*cos(2*pi*F10*t+pi*P10/180); 63 signal11=A11*cos(2*pi*F11*t+pi*P11/180); 64 65 66 % 原始信號+隨機噪音 67 % Y1 = wgn(1,length(signal1),0) + awgn(signal1,10,0); 68 % Y2 = wgn(1,length(signal2),2)+ awgn(signal2,4,5,'linear'); 69 % Y3 = wgn(1,length(signal3),3) + awgn(signal3,4,'measured','linear'); % 70 % Y4 = wgn(1,length(signal4),4) + signal4; 71 72 % 合成信號 73 S = A0+ signal1 + signal2 + signal3 + signal4 + signal5 +signal6+signal7+signal8+signal9; 74 % S = A0+ signal1 + signal2 + signal3 + signal7 + signal10; 75 % S =A0+ signal1 + signal2 + signal3 + signal4 + signal5+ signal6 + signal7 + signal8 + signal9 + signal10; 76 77 78 79 Data = data'; 80 S_data = Data(1:N); 81 82 t2 = 1:1:length(S_data); 83 t1 = 1:1:length(S); 84 85 figure(10); 86 plot(t2, S_data, t1, S); 87 title("頻率為表1中前9組的震動信號"); 88 legend( '靜止-空載的震動信號', '模擬的震動信號'); 89 xlabel('時間'); 90 ylabel('幅值'); 91 92 % 加噪音的合成信號 93 % S = Y1 + Y2 + Y3 + Y4; 94 % S = signal1 + signal2 + signal5 + signal6; 95 96 97 %顯示原始信號 98 figure(1); 99 % plot(t,S1,'b', t, S, 'r'); 100 plot(S); 101 title('混合余弦+正弦信號'); 102 legend('原始的混合信號', '加噪音的信號'); 103 xlabel('時間'); 104 ylabel('幅值'); 105 106 %% 頻率與幅值 107 Y = fft(S, N); 108 Ayy = (abs(Y)); 109 110 % 因為MATLAB中FFT的變換矩陣不是一個酉矩陣(Unitary Matrix)，該陣除以1/sqrt(N)就是個酉矩陣。故經過變換后對信號有放大作用， 111 % 所以要在fft處理后結果除以N/2來，修正此“放大”作用。但是結果在直流的那一點是錯誤的，實際上直流應該除以N修正。 112 F = ([1:N]-1)*Fs/N; % 換成實際的頻率值， N*Fs/2 113 Ayy = Ayy / (N /2); %換算成實際的幅度 114 Ayy(1) = Ayy(1) /2; 115 116 figure(2); 117 plot(F(1:N/2),Ayy(1:N/2)); 118 title('實際幅度-頻率曲線圖'); 119 xlabel('頻率f/Hz'); 120 ylabel('幅度值'); 121 122 %% 頻率與相位 123 % 相位的計算可用函數atan2(b,a)計算。atan2(b,a)是求坐標為(a,b)點的角度值，范圍從-pi到pi。 124 % atan2(500, 148)=x,結果是弧度，換算為角度就是180*(-x)/pi=相位值。 125 Pyy = [1 : N/2]; 126 for i = 1 : N/2 127 Pyy(i) =phase(Y(i)); %計算相位 128 Pyy(i) = Pyy(i) * 180 /pi; %換算為角度 129 end; 130 figure(3); 131 plot(F(1 : N/2), Pyy(1 : N/2)); %顯示相位圖 132 title('相位-頻率曲線圖'); 133 xlabel('頻率f/Hz'); 134 ylabel('初始相位值（角度）');
C++實現：
AudioFFT.h
1 // ================================================================================== 2 // Copyright (c) 2017 HiFi-LoFi 3 // 4 // Permission is hereby granted, free of charge, to any person obtaining a copy 5 // of this software and associated documentation files (the "Software"), to deal 6 // in the Software without restriction, including without limitation the rights 7 // to use, copy, modify, merge, publish, distribute, sublicense, and/or sell 8 // copies of the Software, and to permit persons to whom the Software is furnished 9 // to do so, subject to the following conditions: 10 // 11 // The above copyright notice and this permission notice shall be included in 12 // all copies or substantial portions of the Software. 13 // 14 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 15 // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS 16 // FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR 17 // COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER 18 // IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION 19 // WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. 20 // ================================================================================== 21 22 #ifndef _AUDIOFFT_H 23 #define _AUDIOFFT_H 24 25 26 /** 27 * AudioFFT provides real-to-complex/complex-to-real FFT routines. 28 * 29 * Features: 30 * 31 * - Real-complex FFT and complex-real inverse FFT for power-of-2-sized real data. 32 * 33 * - Uniform interface to different FFT implementations (currently Ooura, FFTW3 and Apple Accelerate). 34 * 35 * - Complex data is handled in "split-complex" format, i.e. there are separate 36 * arrays for the real and imaginary parts which can be useful for SIMD optimizations 37 * (split-complex arrays have to be of length (size/2+1) representing bins from DC 38 * to Nyquist frequency). 39 * 40 * - Output is "ready to use" (all scaling etc. is already handled internally). 41 * 42 * - No allocations/deallocations after the initialization which makes it usable 43 * for real-time audio applications (that's what I wrote it for and using it). 44 * 45 * 46 * How to use it in your project: 47 * 48 * - Add the .h and .cpp file to your project - that's all. 49 * 50 * - To get extra speed, you can link FFTW3 to your project and define 51 * AUDIOFFT_FFTW3 (however, please check whether your project suits the 52 * according license). 53 * 54 * - To get the best speed on Apple platforms, you can link the Apple 55 * Accelerate framework to your project and define 56 * AUDIOFFT_APPLE_ACCELERATE (however, please check whether your 57 * project suits the according license). 58 * 59 * 60 * Remarks: 61 * 62 * - AudioFFT is not intended to be the fastest FFT, but to be a fast-enough 63 * FFT suitable for most audio applications. 64 * 65 * - AudioFFT uses the quite liberal MIT license. 66 * 67 * 68 * Example usage: 69 * @code 70 * #include "AudioFFT.h" 71 * 72 * void Example() 73 * { 74 * const size_t fftSize = 1024; // Needs to be power of 2! 75 * 76 * std::vector<float> input(fftSize, 0.0f); 77 * std::vector<float> re(audiofft::AudioFFT::ComplexSize(fftSize)); 78 * std::vector<float> im(audiofft::AudioFFT::ComplexSize(fftSize)); 79 * std::vector<float> output(fftSize); 80 * 81 * audiofft::AudioFFT fft; 82 * fft.init(1024); 83 * fft.fft(input.data(), re.data(), im.data()); 84 * fft.ifft(output.data(), re.data(), im.data()); 85 * } 86 * @endcode 87 */ 88 89 90 #include <cstddef> 91 #include <memory> 92 93 94 namespace audiofft 95 { 96 97 namespace detail 98 { 99 class AudioFFTImpl; 100 } 101 102 103 // ============================================================= 104 105 106 /** 107 * @class AudioFFT 108 * @brief Performs 1D FFTs 109 */ 110 class AudioFFT 111 { 112 public: 113 /** 114 * @brief Constructor 115 */ 116 AudioFFT(); 117 118 AudioFFT(const AudioFFT&) = delete; 119 AudioFFT& operator=(const AudioFFT&) = delete; 120 121 /** 122 * @brief Destructor 123 */ 124 ~AudioFFT(); 125 126 /** 127 * @brief Initializes the FFT object 128 * @param size Size of the real input (must be power 2) 129 */ 130 void init(size_t size); 131 132 /** 133 * @brief Performs the forward FFT 134 * @param data The real input data (has to be of the length as specified in init()) 135 * @param re The real part of the complex output (has to be of length as returned by ComplexSize()) 136 * @param im The imaginary part of the complex output (has to be of length as returned by ComplexSize()) 137 */ 138 void fft(const float* data, float* re, float* im); 139 140 /** 141 * @brief Performs the inverse FFT 142 * @param data The real output data (has to be of the length as specified in init()) 143 * @param re The real part of the complex input (has to be of length as returned by ComplexSize()) 144 * @param im The imaginary part of the complex input (has to be of length as returned by ComplexSize()) 145 */ 146 void ifft(float* data, const float* re, const float* im); 147 148 /** 149 * @brief Calculates the necessary size of the real/imaginary complex arrays 150 * @param size The size of the real data 151 * @return The size of the real/imaginary complex arrays 152 */ 153 static size_t ComplexSize(size_t size); 154 155 private: 156 std::unique_ptr<detail::AudioFFTImpl> _impl; 157 }; 158 159 160 /** 161 * @deprecated 162 * @brief Let's keep an AudioFFTBase type around for now because it has been here already in the 1st version in order to avoid breaking existing code. 163 */ 164 typedef AudioFFT AudioFFTBase; 165 166 } // End of namespace 167 168 #endif // Header guard
AudioFFT.cpp
// ================================================================================== // Copyright (c) 2017 HiFi-LoFi // // Permission is hereby granted, free of charge, to any person obtaining a copy // of this software and associated documentation files (the "Software"), to deal // in the Software without restriction, including without limitation the rights // to use, copy, modify, merge, publish, distribute, sublicense, and/or sell // copies of the Software, and to permit persons to whom the Software is furnished // to do so, subject to the following conditions: // // The above copyright notice and this permission notice shall be included in // all copies or substantial portions of the Software. // // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS // FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR // COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER // IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION // WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. // ================================================================================== #include "AudioFFT.h" #include <cassert> #include <cmath> #include <cstring> #if defined(AUDIOFFT_INTEL_IPP) #define AUDIOFFT_INTEL_IPP_USED #include <ipp.h> #elif defined(AUDIOFFT_APPLE_ACCELERATE) #define AUDIOFFT_APPLE_ACCELERATE_USED #include <Accelerate/Accelerate.h> #include <vector> #elif defined (AUDIOFFT_FFTW3) #define AUDIOFFT_FFTW3_USED #include <fftw3.h> #else #if !defined(AUDIOFFT_OOURA) #define AUDIOFFT_OOURA #endif #define AUDIOFFT_OOURA_USED #include <vector> #endif namespace audiofft { namespace detail { class AudioFFTImpl { public: AudioFFTImpl() = default; AudioFFTImpl(const AudioFFTImpl&) = delete; AudioFFTImpl& operator=(const AudioFFTImpl&) = delete; virtual ~AudioFFTImpl() = default; virtual void init(size_t size) = 0; virtual void fft(const float* data, float* re, float* im) = 0; virtual void ifft(float* data, const float* re, const float* im) = 0; }; constexpr bool IsPowerOf2(size_t val) { return (val == 1 || (val & (val - 1)) == 0); } template<typename TypeDest, typename TypeSrc> void ConvertBuffer(TypeDest* dest, const TypeSrc* src, size_t len) { for (size_t i = 0; i<len; ++i) { dest[i] = static_cast<TypeDest>(src[i]); } } template<typename TypeDest, typename TypeSrc, typename TypeFactor> void ScaleBuffer(TypeDest* dest, const TypeSrc* src, const TypeFactor factor, size_t len) { for (size_t i = 0; i<len; ++i) { dest[i] = static_cast<TypeDest>(static_cast<TypeFactor>(src[i]) * factor); } } } // End of namespace detail // ================================================================ #ifdef AUDIOFFT_OOURA_USED /** * @internal * @class OouraFFT * @brief FFT implementation based on the great radix-4 routines by Takuya Ooura */ class OouraFFT : public detail::AudioFFTImpl { public: OouraFFT() : detail::AudioFFTImpl(), _size(0), _ip(), _w(), _buffer() { } OouraFFT(const OouraFFT&) = delete; OouraFFT& operator=(const OouraFFT&) = delete; virtual void init(size_t size) override { if (_size != size) { _ip.resize(2 + static_cast<int>(std::sqrt(static_cast<double>(size)))); _w.resize(size / 2); _buffer.resize(size); _size = size; const int size4 = static_cast<int>(_size) / 4; makewt(size4, _ip.data(), _w.data()); makect(size4, _ip.data(), _w.data() + size4); } } virtual void fft(const float* data, float* re, float* im) override { // Convert into the format as required by the Ooura FFT detail::ConvertBuffer(_buffer.data(), data, _size); rdft(static_cast<int>(_size), +1, _buffer.data(), _ip.data(), _w.data()); // Convert back to split-complex { double* b = _buffer.data(); double* bEnd = b + _size; float *r = re; float *i = im; while (b != bEnd) { *(r++) = static_cast<float>(*(b++)); *(i++) = static_cast<float>(-(*(b++))); } } const size_t size2 = _size / 2; re[size2] = -im[0]; im[0] = 0.0; im[size2] = 0.0; } virtual void ifft(float* data, const float* re, const float* im) override { // Convert into the format as required by the Ooura FFT { double* b = _buffer.data(); double* bEnd = b + _size; const float *r = re; const float *i = im; while (b != bEnd) { *(b++) = static_cast<double>(*(r++)); *(b++) = -static_cast<double>(*(i++)); } _buffer[1] = re[_size / 2]; } rdft(static_cast<int>(_size), -1, _buffer.data(), _ip.data(), _w.data()); // Convert back to split-complex detail::ScaleBuffer(data, _buffer.data(), 2.0 / static_cast<double>(_size), _size); } private: size_t _size; std::vector<int> _ip; std::vector<double> _w; std::vector<double> _buffer; void rdft(int n, int isgn, double *a, int *ip, double *w) { int nw = ip[0]; int nc = ip[1]; if (isgn >= 0) { if (n > 4) { bitrv2(n, ip + 2, a); cftfsub(n, a, w); rftfsub(n, a, nc, w + nw); } else if (n == 4) { cftfsub(n, a, w); } double 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 + 2, a); cftbsub(n, a, w); } else if (n == 4) { cftfsub(n, a, w); } } } /* -------- initializing routines -------- */ 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 = atan(1.0) / nwh; w[0] = 1; w[1] = 0; w[nwh] = cos(delta * nwh); w[nwh + 1] = w[nwh]; if (nwh > 2) { for (j = 2; j < nwh; j += 2) { x = cos(delta * j); y = sin(delta * j); w[j] = x; w[j + 1] = y; w[nw - j] = y; w[nw - j + 1] = x; } bitrv2(nw, ip + 2, w); } } } void makect(int nc, int *ip, double *c) { int j, nch; double delta; ip[1] = nc; if (nc > 1) { nch = nc >> 1; delta = atan(1.0) / nch; c[0] = cos(delta * nch); c[nch] = 0.5 * c[0]; for (j = 1; j < nch; j++) { c[j] = 0.5 * cos(delta * j); c[nc - j] = 0.5 * sin(delta * j); } } } /* -------- child routines -------- */ void bitrv2(int n, int *ip, double *a) { int j, j1, k, k1, l, m, m2; double xr, xi, yr, yi; ip[0] = 0; l = n; m = 1; while ((m << 3) < l) { l >>= 1; for (j = 0; j < m; j++) { ip[m + j] = ip[j] + 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]; k1 = 2 * k + ip[j]; 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]; 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]; k1 = 2 * k + ip[j]; 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; } } } } 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; } } } 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; } } } 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; } } 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; } } } void rftfsub(int n, double *a, int nc, double *c) { 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[nc - kk]; wki = c[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; } } void rftbsub(int n, double *a, int nc, double *c) { 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[nc - kk]; wki = c[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]; } }; /** * @internal * @brief Concrete FFT implementation */ typedef OouraFFT AudioFFTImplementation; #endif // AUDIOFFT_OOURA_USED // ================================================================ #ifdef AUDIOFFT_INTEL_IPP_USED /** * @internal * @class IntelIppFFT * @brief FFT implementation using the Intel Integrated Performance Primitives */ class IntelIppFFT : public detail::AudioFFTImpl { public: IntelIppFFT() : detail::AudioFFTImpl(), _size(0), _operationalBufferSize(0), _powerOf2(0), _fftSpec(nullptr), _fftSpecBuf(0), _fftWorkBuf(0), _operationalBuffer(nullptr) { ippInit(); } IntelIppFFT(const IntelIppFFT&) = delete; IntelIppFFT& operator=(const IntelIppFFT&) = delete; virtual ~IntelIppFFT() { init(0); } virtual void init(size_t size) override { if (_fftSpec) { if (_fftWorkBuf) ippFree(_fftWorkBuf); if (_fftSpecBuf) ippFree(_fftSpecBuf); ippFree(_operationalBuffer); _size = 0; _operationalBufferSize = 0; _powerOf2 = 0; _fftSpec = 0; } if (size > 0) { _size = size; _operationalBufferSize = _size + 2; _powerOf2 = (int)(log((double)_size) / log(2.0)); // Query to get buffer sizes int sizeFFTSpec, sizeFFTInitBuf, sizeFFTWorkBuf; ippsFFTGetSize_R_32f( _powerOf2, IPP_FFT_NODIV_BY_ANY, ippAlgHintAccurate, &sizeFFTSpec, &sizeFFTInitBuf, &sizeFFTWorkBuf ); Ipp8u* fftInitBuf; // init buffers _fftSpecBuf = ippsMalloc_8u(sizeFFTSpec); _fftWorkBuf = ippsMalloc_8u(sizeFFTWorkBuf); fftInitBuf = ippsMalloc_8u(sizeFFTInitBuf); // Initialize FFT ippsFFTInit_R_32f( &_fftSpec, _powerOf2, IPP_FFT_NODIV_BY_ANY, ippAlgHintAccurate, _fftSpecBuf, fftInitBuf ); if (fftInitBuf) ippFree(fftInitBuf); // init operational buffer _operationalBuffer = ippsMalloc_32f( _operationalBufferSize ); } } virtual void fft(const float* data, float* re, float* im) override { size_t complexNumbersCount = _operationalBufferSize / 2; ippsFFTFwd_RToCCS_32f( data, _operationalBuffer, _fftSpec, _fftWorkBuf ); // no need to scale size_t complexCounter = 0; for (int i = 0; i < complexNumbersCount; ++i) { re[i] = _operationalBuffer[complexCounter++]; im[i] = _operationalBuffer[complexCounter++]; } } virtual void ifft(float* data, const float* re, const float* im) override { size_t complexNumbersCount = _operationalBufferSize / 2; size_t complexCounter = 0; for (int i = 0; i < complexNumbersCount; ++i) { _operationalBuffer[complexCounter++] = re[i]; _operationalBuffer[complexCounter++] = im[i]; } ippsFFTInv_CCSToR_32f( _operationalBuffer, data, _fftSpec, _fftWorkBuf ); // scaling const float factor = 1.0f / static_cast<float>(_size); ippsMulC_32f_I(factor, data, _size); } private: size_t _size; size_t _operationalBufferSize; size_t _powerOf2; IppsFFTSpec_R_32f* _fftSpec; Ipp8u* _fftSpecBuf; Ipp8u* _fftWorkBuf; Ipp32f* _operationalBuffer; }; /** * @internal * @brief Concrete FFT implementation */ typedef IntelIppFFT AudioFFTImplementation; #endif // AUDIOFFT_INTEL_IPP_USED // ================================================================ #ifdef AUDIOFFT_APPLE_ACCELERATE_USED /** * @internal * @class AppleAccelerateFFT * @brief FFT implementation using the Apple Accelerate framework internally */ class AppleAccelerateFFT : public detail::AudioFFTImpl { public: AppleAccelerateFFT() : detail::AudioFFTImpl(), _size(0), _powerOf2(0), _fftSetup(0), _re(), _im() { } AppleAccelerateFFT(const AppleAccelerateFFT&) = delete; AppleAccelerateFFT& operator=(const AppleAccelerateFFT&) = delete; virtual ~AppleAccelerateFFT() { init(0); } virtual void init(size_t size) override { if (_fftSetup) { vDSP_destroy_fftsetup(_fftSetup); _size = 0; _powerOf2 = 0; _fftSetup = 0; _re.clear(); _im.clear(); } if (size > 0) { _size = size; _powerOf2 = 0; while ((1 << _powerOf2) < _size) { ++_powerOf2; } _fftSetup = vDSP_create_fftsetup(_powerOf2, FFT_RADIX2); _re.resize(_size / 2); _im.resize(_size / 2); } } virtual void fft(const float* data, float* re, float* im) override { const size_t size2 = _size / 2; DSPSplitComplex splitComplex; splitComplex.realp = re; splitComplex.imagp = im; vDSP_ctoz(reinterpret_cast<const COMPLEX*>(data), 2, &splitComplex, 1, size2); vDSP_fft_zrip(_fftSetup, &splitComplex, 1, _powerOf2, FFT_FORWARD); const float factor = 0.5f; vDSP_vsmul(re, 1, &factor, re, 1, size2); vDSP_vsmul(im, 1, &factor, im, 1, size2); re[size2] = im[0]; im[0] = 0.0f; im[size2] = 0.0f; } virtual void ifft(float* data, const float* re, const float* im) override { const size_t size2 = _size / 2; ::memcpy(_re.data(), re, size2 * sizeof(float)); ::memcpy(_im.data(), im, size2 * sizeof(float)); _im[0] = re[size2]; DSPSplitComplex splitComplex; splitComplex.realp = _re.data(); splitComplex.imagp = _im.data(); vDSP_fft_zrip(_fftSetup, &splitComplex, 1, _powerOf2, FFT_INVERSE); vDSP_ztoc(&splitComplex, 1, reinterpret_cast<COMPLEX*>(data), 2, size2); const float factor = 1.0f / static_cast<float>(_size); vDSP_vsmul(data, 1, &factor, data, 1, _size); } private: size_t _size; size_t _powerOf2; FFTSetup _fftSetup; std::vector<float> _re; std::vector<float> _im; }; /** * @internal * @brief Concrete FFT implementation */ typedef AppleAccelerateFFT AudioFFTImplementation; #endif // AUDIOFFT_APPLE_ACCELERATE_USED // ================================================================ #ifdef AUDIOFFT_FFTW3_USED /** * @internal * @class FFTW3FFT * @brief FFT implementation using FFTW3 internally (see fftw.org) */ class FFTW3FFT : public detail::AudioFFTImpl { public: FFTW3FFT() : detail::AudioFFTImpl(), _size(0), _complexSize(0), _planForward(0), _planBackward(0), _data(0), _re(0), _im(0) { } FFTW3FFT(const FFTW3FFT&) = delete; FFTW3FFT& operator=(const FFTW3FFT&) = delete; virtual ~FFTW3FFT() { init(0); } virtual void init(size_t size) override { if (_size != size) { if (_size > 0) { fftwf_destroy_plan(_planForward); fftwf_destroy_plan(_planBackward); _planForward = 0; _planBackward = 0; _size = 0; _complexSize = 0; if (_data) { fftwf_free(_data); _data = 0; } if (_re) { fftwf_free(_re); _re = 0; } if (_im) { fftwf_free(_im); _im = 0; } } if (size > 0) { _size = size; _complexSize = AudioFFT::ComplexSize(_size); const size_t complexSize = AudioFFT::ComplexSize(_size); _data = reinterpret_cast<float*>(fftwf_malloc(_size * sizeof(float))); _re = reinterpret_cast<float*>(fftwf_malloc(complexSize * sizeof(float))); _im = reinterpret_cast<float*>(fftwf_malloc(complexSize * sizeof(float))); fftw_iodim dim; dim.n = static_cast<int>(size); dim.is = 1; dim.os = 1; _planForward = fftwf_plan_guru_split_dft_r2c(1, &dim, 0, 0, _data, _re, _im, FFTW_MEASURE); _planBackward = fftwf_plan_guru_split_dft_c2r(1, &dim, 0, 0, _re, _im, _data, FFTW_MEASURE); } } } virtual void fft(const float* data, float* re, float* im) override { ::memcpy(_data, data, _size * sizeof(float)); fftwf_execute_split_dft_r2c(_planForward, _data, _re, _im); ::memcpy(re, _re, _complexSize * sizeof(float)); ::memcpy(im, _im, _complexSize * sizeof(float)); } virtual void ifft(float* data, const float* re, const float* im) override { ::memcpy(_re, re, _complexSize * sizeof(float)); ::memcpy(_im, im, _complexSize * sizeof(float)); fftwf_execute_split_dft_c2r(_planBackward, _re, _im, _data); detail::ScaleBuffer(data, _data, 1.0f / static_cast<float>(_size), _size); } private: size_t _size; size_t _complexSize; fftwf_plan _planForward; fftwf_plan _planBackward; float* _data; float* _re; float* _im; }; /** * @internal * @brief Concrete FFT implementation */ typedef FFTW3FFT AudioFFTImplementation; #endif // AUDIOFFT_FFTW3_USED // ============================================================= AudioFFT::AudioFFT() : _impl(new AudioFFTImplementation()) { } AudioFFT::~AudioFFT() { } void AudioFFT::init(size_t size) { assert(detail::IsPowerOf2(size)); _impl->init(size); } void AudioFFT::fft(const float* data, float* re, float* im) { _impl->fft(data, re, im); } void AudioFFT::ifft(float* data, const float* re, const float* im) { _impl->ifft(data, re, im); } size_t AudioFFT::ComplexSize(size_t size) { return (size / 2) + 1; } } // End of namespace
FrequencyAnalyzer.h
1 #pragma once 2 3 #include <vector> 4 5 6 typedef struct { 7 float frequency, amplitude, phase; 8 int sin_cos_flag = 1; // 0: 正弦 1:余弦 9 }FrequencyInfo; 10 11 typedef struct { 12 13 //自行定義參數 14 float sampleFrequency; // 采樣頻率 15 int dataLength; // 頻域分析處理的數據長度 16 17 }Cfg; 18 19 class CFrequencyAnalyzer { 20 21 22 public: 23 CFrequencyAnalyzer() {}; 24 virtual ~CFrequencyAnalyzer() {}; 25 26 /* 27 功能：分析信號頻域信息(賦值/頻率/相位) 28 形參： 29 signal:稱重傳感器信號 30 返回值：信號的頻域信息 31 */ 32 virtual std::vector<FrequencyInfo> calcSignalFrequency(const std::vector<int>& signal) = 0; 33 34 // 合成各子信號 35 /* 36 功能：合成各子信號(賦值/頻率/相位) 37 形參： 38 list_signal: 各子信號參數 39 begin:起始位置 40 Fs: 采樣頻率 41 size:采樣點數 42 返回值：合成后的信號 43 */ 44 virtual std::vector<int> calcSynthesizedSignal(const std::vector<FrequencyInfo>& list_signal, const int begin, const int Fs, const int size) = 0; 45 46 // 功能：檢查參數合法性，若合法則接收保存 47 virtual int setCfg(const Cfg& cfg) = 0; 48 49 // 功能：返回參數 50 const Cfg& getCfg() { return _cfg; } 51 52 protected: 53 Cfg _cfg; 54 };
FrequencyAnalyzerDFT.h
#ifndef _FREQEUENCY_ANLAYZER_DFT_H_DEF_ #define _FREQEUENCY_ANLAYZER_DFT_H_DEF_ #include "FrequencyAnalyzer.h" #include "AudioFFT.h" #define pi 3.1415926 class FrequencyAnalyzerDFT : public CFrequencyAnalyzer { public: /******************************************************************************* 功能：構造函數 *******************************************************************************/ FrequencyAnalyzerDFT(); /******************************************************************************* 功能：析構函數 *******************************************************************************/ virtual ~FrequencyAnalyzerDFT(); /* 功能：分析信號頻域信息(賦值/頻率/相位) 形參： signal:稱重傳感器信號 返回值：信號的頻域信息 */ virtual std::vector<FrequencyInfo> calcSignalFrequency(const std::vector<int>& signal); // 合成各子信號 /* 功能：合成各子信號(賦值/頻率/相位) 形參： list_signal: 各子信號參數 begin:起始位置 Fs: 采樣頻率 size:采樣點數 返回值：合成后的信號 */ virtual std::vector<int> calcSynthesizedSignal(const std::vector<FrequencyInfo>& list_signal, const int begin, const int Fs, const int size); // 功能：檢查參數合法性，若合法則接收保存 virtual int setCfg(const Cfg& cfg); // 功能：返回參數 const Cfg& getCfg() { return _cfg; } private: Cfg _cfg; audiofft::AudioFFT fft; int m_dataLength; // 數據長度 float m_sampleFrequency; // 采樣頻率 std::vector<float> m_dataInput; // 輸入數據 std::vector<float> m_fft_real; // fft后的實部 std::vector<float> m_fft_imaginary; // fft后的虛部 }; #endif //_FREQEUENCY_ANLAYZER_DFT_H_DEF_
FrequencyAnalyzerDFT.cpp
#include "FrequencyAnalyzerDFT.h" FrequencyAnalyzerDFT::FrequencyAnalyzerDFT() { } FrequencyAnalyzerDFT::~FrequencyAnalyzerDFT() { } std::vector<FrequencyInfo> FrequencyAnalyzerDFT::calcSignalFrequency(const std::vector<int>& signal) { int length = signal.size(); std::vector<FrequencyInfo> frequencyInfo; if (length < _cfg.dataLength) { return frequencyInfo; } for (int i = 0; i < _cfg.dataLength; i++) { m_dataInput[i] = (float)signal[i]; } fft.fft(m_dataInput.data(), m_fft_real.data(), m_fft_imaginary.data()); int fftHalfSize = _cfg.dataLength / 2; float fs_per = _cfg.sampleFrequency / _cfg.dataLength; float ampTemp = _cfg.dataLength / 2.0; float Amplitude = 0.0, Frequency = 0.0, Phase = 0.0; FrequencyInfo f_info; // 計算頻率與幅值 for (int m = 0; m < fftHalfSize; m++){ // 頻率 Frequency= m * fs_per; // 幅值 Amplitude = abs(sqrt(pow(m_fft_real[m], 2) + pow(m_fft_imaginary[m], 2))); Amplitude = Amplitude / ampTemp; // 頻率為0時，幅值減半 if (m == 0) { Amplitude = Amplitude / 2.0; } // 初始相位 Phase = atan2f(m_fft_imaginary[m], m_fft_real[m]); Phase = Phase * 180 / pi; f_info.amplitude = Amplitude; f_info.phase = Phase; f_info.frequency = Frequency; // 余弦函數的標志 f_info.sin_cos_flag = 1; frequencyInfo.push_back(f_info); } return frequencyInfo; } std::vector<int> FrequencyAnalyzerDFT::calcSynthesizedSignal(const std::vector<FrequencyInfo>& list_signal, const int begin, const int Fs, const int size) { std::vector<int> synData; if (size < 0 || list_signal.size()<0) return synData; int startIndex = begin; int endIndex = begin + size; std::vector<std::vector<float>> allData; for (int t = 0; t < endIndex; t++) { std::vector<float> singleData; float temp = 0.0; for (int i = 0; i < list_signal.size(); i++) { float A = list_signal[i].amplitude; float F = list_signal[i].frequency; float P = list_signal[i].phase; int flag = list_signal[i].sin_cos_flag; // printf("%f %f %f \n", F, A, P); if (flag == 1) { temp += A*cos(2.0 * pi*F / Fs*t + P*pi / 180.0); } else { temp += A*sin(2.0 * pi*F / Fs*t + P*pi / 180.0); } } synData.push_back(int(temp)); } return synData; } int FrequencyAnalyzerDFT::setCfg(const Cfg& cfg) { int errorCode = 0; if (cfg.dataLength < 0) { errorCode = -1; return errorCode; } if (cfg.sampleFrequency < 0){ errorCode = -1; return errorCode; } _cfg.dataLength = cfg.dataLength; _cfg.sampleFrequency = cfg.sampleFrequency; m_dataLength = _cfg.dataLength; m_sampleFrequency = _cfg.sampleFrequency; printf("%d, %f \n", m_dataLength, m_sampleFrequency); fft.init(m_dataLength); m_dataInput.assign(_cfg.dataLength, 0.0f); m_fft_real.assign(audiofft::AudioFFT::ComplexSize(m_dataLength), 0.0f); m_fft_imaginary.assign(audiofft::AudioFFT::ComplexSize(m_dataLength), 0.0f); return errorCode; }
main.cpp
1 #include <iostream> 2 #include <fstream> 3 #include <iomanip> 4 #include <Windows.h> 5 6 #include "AudioFFT.h" 7 #include "FrequencyAnalyzerDFT.h" 8 #include "FrequencyAnalyzer.h" 9 10 #include <string> 11 #include <vector> 12 #include <complex> 13 14 using namespace audiofft; 15 using namespace std; 16 17 #define pi 3.1415926 18 19 void LoadData(std::string &path, int &rows, int &cols, std::vector<std::vector<complex<double>>> &Data) 20 { 21 22 //ifstream myfile("E:\\新北洋_文檔\\動態稱\\傅里葉變換的耗時驗證\\DATA\\16_自-20kg-400\\demo.txt"); 23 ifstream myfile(path); 24 25 if (!myfile.is_open()) 26 { 27 cout << "can not open this file" << endl; 28 return; 29 } 30 31 //std::vector<std::vector<complex<double>>> Data; 32 std::vector<complex<double>> data; 33 34 for (int i = 0; i < cols; i++) 35 { 36 data.clear(); 37 38 for (int j = 0; j < rows; j++) 39 { 40 complex<double> dataTemp; 41 42 myfile >> dataTemp; 43 data.push_back(dataTemp); 44 } 45 46 Data.push_back(data); 47 } 48 49 myfile.close(); 50 } 51 52 void AudioFFTTest(std::vector<std::vector<complex<double>>> &Data, double &Fs) 53 { 54 55 if (Data.size() < 0) 56 return; 57 58 double fs = Fs; 59 double Wc = 20; // 截止頻率 60 61 int pre = 18; // 保存的數據精度 62 63 ofstream outFile1("fft_re.txt"); // 實部 64 65 66 ofstream outFile2("fft_im.txt"); // 虛部 67 ofstream outFile3("out.txt"); // data -> FFT --> ifft --> dataOut 68 69 ofstream outFile4("Wc.txt"); // 頻率 70 ofstream outFile5("A.txt"); // 幅值 71 ofstream outFile6("Angle.txt"); // 72 73 outFile1.precision(pre); 74 outFile2.precision(pre); 75 outFile3.precision(pre); 76 outFile4.precision(pre); 77 outFile5.precision(pre); 78 outFile6.precision(pre); 79 80 for (int i = 0; i < Data.size(); i++) 81 { 82 const size_t fftSize = Data[i].size(); 83 std::vector<float> input(fftSize, 0.0f); 84 std::vector<float> re(audiofft::AudioFFT::ComplexSize(fftSize)); 85 std::vector<float> im(audiofft::AudioFFT::ComplexSize(fftSize)); 86 std::vector<float> output(fftSize); 87 88 for (int j = 0; j < fftSize; j++) 89 { 90 input[j] = (float)Data[i][j].real(); 91 } 92 93 audiofft::AudioFFT fft; 94 fft.init(fftSize); 95 96 // fft 97 fft.fft(input.data(), re.data(), im.data()); 98 99 int fftHalfSize = fftSize / 2; 100 101 double wc_temp = 0, Amp_temp = 0; 102 double fs_per = fs / fftSize; 103 104 std::vector<double> Amplitude, Frequency; 105 106 double ampTemp = fftSize / 2.0; 107 108 // 計算頻率與幅值 109 for (int m = 0; m <fftHalfSize; m++) 110 { 111 wc_temp = (double)(m*fs_per); 112 113 Amp_temp = abs(sqrt(pow(re[m], 2) + pow(im[m], 2))); 114 Amp_temp = Amp_temp / ampTemp; 115 116 Frequency.push_back(wc_temp); 117 Amplitude.push_back(Amp_temp); 118 } 119 120 Amplitude[0] = Amplitude[0] / 2; 121 122 std::vector<double> angleVec; 123 double angle = 0.0; 124 for (int j = 0; j < fftHalfSize; j++) 125 { 126 angle = atan2f( im[j], re[j]); 127 angle = angle * 180 / pi; 128 angleVec.push_back(angle); 129 } 130 131 for (int j = 0; j < fftHalfSize; j++) 132 { 133 outFile4 << Frequency[j] << endl; 134 outFile5 << Amplitude[j] << endl; 135 outFile6 << angleVec[j] << endl; 136 } 137 138 // ifft 139 fft.ifft(output.data(), re.data(), im.data()); 140 141 for (int j = 0; j < fftSize; j++) 142 { 143 printf("%10f\n ", re[j]); 144 outFile1 << re[j]<< endl; 145 outFile2 << im[j] << endl; 146 outFile3 << output[j] << endl; 147 } 148 149 } 150 151 outFile1.close(); 152 outFile2.close(); 153 outFile3.close(); 154 outFile4.close(); 155 outFile5.close(); 156 outFile6.close(); 157 } 158 159 void singleData(std::vector<std::vector<complex<double>>> Data, std::vector<int> &inputData) 160 { 161 if (Data.size() < 0) 162 return; 163 164 for (int i = 0; i < 1; i++) 165 { 166 for (int j = 0; j < Data[0].size(); j++) 167 { 168 inputData.push_back((int)Data[i][j].real()); 169 } 170 } 171 172 } 173 174 void process() 175 { 176 177 std::string path = "E:\\Project\\Matlab_Pro\\DynamicBalance\\Data.txt"; 178 179 int rows = 256, cols = 1; 180 std::vector<std::vector<complex<double>>> Data; 181 182 // 加載數據 183 LoadData(path, rows, cols, Data); 184 185 /*double Fs = 1200; 186 AudioFFTTest(Data, Fs); 187 188 return;*/ 189 190 std::vector<int> inputData; 191 singleData(Data, inputData); 192 193 int dataLength = 256; 194 float sampleFrequency = 1200; 195 196 FrequencyAnalyzerDFT* freqAnalyDft; 197 freqAnalyDft = new FrequencyAnalyzerDFT(); 198 199 Cfg cfg = { 1200, 256 }; 200 //cfg.dataLength = 256; 201 //cfg.sampleFrequency = 1200; 202 203 freqAnalyDft->setCfg(cfg); 204 205 std::vector<FrequencyInfo> frequencyInfo; 206 frequencyInfo = freqAnalyDft->calcSignalFrequency(inputData); 207 208 for (int i = 0; i < frequencyInfo.size(); i++) 209 { 210 printf("%f %f %f %d \n", frequencyInfo[i].frequency, frequencyInfo[i].amplitude, 211 frequencyInfo[i].phase, frequencyInfo[i].sin_cos_flag); 212 } 213 214 std::vector<int> sysData; 215 sysData = freqAnalyDft->calcSynthesizedSignal(frequencyInfo, 0, 1200, 256); 216 217 for (int i = 0; i < sysData.size(); i++) 218 { 219 printf("%d \n", sysData[i]); 220 } 221 222 } 223 224 int main() 225 { 226 227 process(); 228 229 system("pause"); 230 231 return 0; 232 }