NOVINKA - Online rekvalifikační kurz Python programátor. Oblíbená a studenty ověřená rekvalifikace - nyní i online.
Hledáme nové posily do ITnetwork týmu. Podívej se na volné pozice a přidej se do nejagilnější firmy na trhu - Více informací.

Diskuze: Výpočet FFT a drobná úprava programu

V předchozím kvízu, Online test znalostí C++, jsme si ověřili nabyté zkušenosti z kurzu.

Aktivity
Avatar
Caster
Člen
Avatar
Caster:23.1.2023 21:36

Pro 32-bitový MCU bych rád upravil C++ program pro výpočet FFT, který jsem našel na geeksforgeeks. Header stdc++.h jsem zkopíroval z odkazu zde, přidal ho do adresáře C:\Program Files\Microsoft Visual Studio\2022\Com­munity\VC\Auxi­liary\VS\inclu­de na svém notebooku a restartoval VS2022. Ještě jsem musel nadefinovat proměnnou M_PI

Zkusil jsem: Program jsem úspěšně přeložil a spustil, funguje.

#include <iostream>

#include <vector>
#include <stdc++.h>
using namespace std;

#define M_PI    3.141592653589793

// For storing complex values of nth roots
// of unity we use complex<double>
typedef complex<double> cd;

// Recursive function of FFT
vector<cd> fft(vector<cd>& a)
{
        int n = a.size();

        // if input contains just one element
        if (n == 1)
                return vector<cd>(1, a[0]);

        // For storing n complex nth roots of unity
        vector<cd> w(n);
        for (int i = 0; i < n; i++) {
                double alpha = -2 * M_PI * i / n;
                w[i] = cd(cos(alpha), sin(alpha));
        }

        vector<cd> A0(n / 2), A1(n / 2);
        for (int i = 0; i < n / 2; i++) {

                // even indexed coefficients
                A0[i] = a[i * 2];

                // odd indexed coefficients
                A1[i] = a[i * 2 + 1];
        }

        // Recursive call for even indexed coefficients
        vector<cd> y0 = fft(A0);

        // Recursive call for odd indexed coefficients
        vector<cd> y1 = fft(A1);

        // for storing values of y0, y1, y2, ..., yn-1.
        vector<cd> y(n);

        for (int k = 0; k < n / 2; k++) {
                y[k] = y0[k] + w[k] * y1[k];
                y[k + n / 2] = y0[k] - w[k] * y1[k];
        }
        return y;
}

// Driver code
int main()
{
        vector<cd> a{ 1, 2, 3, 4 };
        vector<cd> b = fft(a);
        for (int i = 0; i < 4; i++)
                cout << b[i] << endl;

        return 0;
}

Chci docílit: Potřebuji poradit, jak upravit výpis, aby našel pouze maximální hodnotu (frekvenci) sqrt (real * real + imag * imag) z vypočtených hodnot s tím, že se začne hledat od 1ní hodnoty, 0tá hodnota je pouze stejnosměrný ofset, jak je vysvětleno v odkazu zde.

V dalším kroku bych chtěl program upravit, aby šel přeložit pro 32-bitové MCU (ATSAMD21J18) na embeded systému.

Pracuji na projektu dálkového měření objemu kapaliny v plastovém septiku 5m3, poklepem solenoidu na vnější stranu septiku, zachycení odezvy pomocí piezoelektrického snímače, výpočtem rezonanční frekvence pomocí FFT což pomocí změřené převodní funkce frekvence na objem umožní zjistit kolik je kapaliny v litrech v nádrži. Výslednou hodnotu pak pomocí sítě LoRaWAN odešlu na mobil.

 
Odpovědět
23.1.2023 21:36
Avatar
DarkCoder
Člen
Avatar
Odpovídá na Caster
DarkCoder:23.1.2023 22:02 
int main()
{
        vector<cd> a{ 1, 2, 3, 4 };
        vector<cd> b = fft(a);
        double maxValue = 0;
        int maxIndex = 0;
        for (int i = 0; i < b.size(); i++)
        {
                double currentValue = sqrt(b[i].real() * b[i].real() + b[i].imag() * b[i].imag());
                if(currentValue > maxValue) {
                   maxValue = currentValue;
                   maxIndex = i;
                }
        }
        cout << "Maximum value (frequency) is: " << maxValue << " found at index: " << maxIndex << endl;
        return 0;
}
Nahoru Odpovědět
23.1.2023 22:02
"I ta nejlepší poučka postrádá na významu, není-li patřičně předána." - DarkCoder
Avatar
Caster
Člen
Avatar
Odpovídá na DarkCoder
Caster:23.1.2023 22:16

Díky za bleskovou reakci. Zkusil jsem to a program vypíše:
Maximum value (frequency) is: 10 found at index: 0

Výsledné hodnoty jsou bez hledání maxima:
(10,0)
(-2,2)
(-2,0)
(-2,-2)

Jak jsem psal, 0tá hodnota je jen stejnosměrný ofset, který se nepočítá. Výsledem by tedy měl být správně asi:
Maximum value (frequency) is: 2,83 found at index: 1

 
Nahoru Odpovědět
23.1.2023 22:16
Avatar
DarkCoder
Člen
Avatar
Odpovídá na Caster
DarkCoder:23.1.2023 22:19

Pak oprav řádek s for cyklem na:

for (int i = 1; i < b.size(); i++)
Nahoru Odpovědět
23.1.2023 22:19
"I ta nejlepší poučka postrádá na významu, není-li patřičně předána." - DarkCoder
Avatar
Caster
Člen
Avatar
Caster:23.1.2023 22:59

Už to počítá správně, díky.
Maximum value (frequency) is: 2.82843 found at index: 1

 
Nahoru Odpovědět
23.1.2023 22:59
Avatar
Caster
Člen
Avatar
Caster:31.1.2023 2:14

Program jsem ještě upravil a dolaďuji ho tak, aby mi přesněji spočítal frekvenci vstupního signálu (4,82 kHz a ne nyní 5 kHz). Po spuštění vypíše frekvenční charakteristiku signálu s vyzdvihnutím hlavní frekvence a potlačením nepodtatných hodnot.

Na třetím řádku výstupu (printf) je max. hodnota signálu a jeho frekvence:
153.005859 5000.000244

Program je níže, pro lepší rozlišení jsem do něj přidal funkci Hanning (použil jsem její kód z MATLAB viz diskuze), která zvýrazní frekvenci signálu a potlačí ostatní (viz popis zde). Potřebuji poradit, jak ji zavolat pro vynásobení hodnot pole xn[] před vlastním výpočtem FFT po úpravě dat pole pod:

i = len;
while (i-- > 0) {
        x[i] = (xn[i] - mean);
}

V MATLAB kódu vypadá její použití takto:
FTs = fft((s-mean(s)).*hann(L),NFFT)/L;
".*" znamená, že funkce hann vynásobí všechny prvky pole s, L je délka pole 1200
NFFT je 2048 (nejbližší větší mocnina násobku 2 délky pole 1200) tj. 211 Ještě přesně nevím, co vlastně počítá NFFT)/L = 2048/1200 = 1,71 Asi to nějak souvisí s upřesněním časové osy signálu.

Více podrobnějších informací je v diskuzi How to interpret results of FFT/DFT?
Vstupní tabulka obsahuje dvě pole. První je jen číslo prvku a druhé je hodnota napětí generované sinusovky o kmitočtu 4,82 kHz (to by měl program v C++ přesně spočítat viz výpočet v MATLABU). V diskuzi jsem později upřesnil, že první pole je vlastně čas hodnot napětí sinusovky s pevným intervalem 0,5 µs (viz také časová základna signálu na osciloskopu 50 µs/čtverec tj. 12x 50 µs = 1200 x 0,5 µs). Cílem mého snažení je připravit, jak budu na 32-bitovém MCU ATSAMD21J18 počítat výslednou frekvenci septiku 5m3 (tj. objem kapliny v něm) po jeho "vybuzení" ťuknutím do vnější strany solenoidem a sejmutím odezvy piezo snímačem, který už je na cestě.

Celý program:

#include <iostream>

#include <math.h>
#include <stdio.h>

#define PI      3.141592653589793

const int len = 1200;

float xn[len] = { 0.284, 0.292, 0.288, 0.276, 0.280, 0.276, 0.276, 0.268, 0.272, 0.264, 0.276, 0.260, 0.264, 0.256, 0.264, 0.256, 0.252, 0.264, 0.260, 0.248, 0.252, 0.248, 0.256, 0.244, 0.248, 0.244, 0.244, 0.248, 0.240, 0.248, 0.244, 0.236, 0.244, 0.236, 0.236, 0.244, 0.240, 0.236, 0.236, 0.244, 0.236, 0.244, 0.244, 0.236, 0.240, 0.240, 0.244, 0.236, 0.236, 0.248, 0.244, 0.236, 0.248, 0.244, 0.248, 0.240, 0.244, 0.248, 0.244, 0.256, 0.244, 0.256, 0.256, 0.248, 0.252, 0.260, 0.260, 0.252, 0.256, 0.264, 0.268, 0.260, 0.264, 0.276, 0.276, 0.268, 0.268, 0.284, 0.280, 0.288, 0.276, 0.288, 0.288, 0.296, 0.288, 0.300, 0.304, 0.296, 0.304, 0.296, 0.308, 0.316, 0.316, 0.308, 0.316, 0.324, 0.328, 0.320, 0.324, 0.340, 0.332, 0.340, 0.336, 0.352, 0.352, 0.344, 0.344, 0.364, 0.364, 0.356, 0.360, 0.376, 0.372, 0.380, 0.376, 0.392, 0.388, 0.392, 0.388, 0.408, 0.408, 0.404, 0.404, 0.424, 0.420, 0.424, 0.420, 0.432, 0.432, 0.440, 0.432, 0.452, 0.456, 0.448, 0.448, 0.468, 0.464, 0.476, 0.464, 0.480, 0.488, 0.480, 0.480, 0.492, 0.496, 0.504, 0.496, 0.508, 0.508, 0.520, 0.520, 0.512, 0.524, 0.536, 0.532, 0.540, 0.532, 0.552, 0.552, 0.544, 0.548, 0.572, 0.560, 0.568, 0.564, 0.584, 0.584, 0.576, 0.580, 0.600, 0.592, 0.600, 0.588, 0.612, 0.608, 0.616, 0.608, 0.628, 0.628, 0.624, 0.616, 0.640, 0.640, 0.632, 0.636, 0.656, 0.656, 0.648, 0.648, 0.668, 0.668, 0.660, 0.660, 0.676, 0.680, 0.672, 0.668, 0.688, 0.684, 0.692, 0.680, 0.700, 0.692, 0.704, 0.704, 0.692, 0.704, 0.712, 0.708, 0.704, 0.712, 0.716, 0.720, 0.712, 0.716, 0.732, 0.720, 0.732, 0.720, 0.740, 0.736, 0.728, 0.732, 0.744, 0.740, 0.732, 0.736, 0.748, 0.740, 0.748, 0.740, 0.752, 0.752, 0.744, 0.744, 0.756, 0.748, 0.756, 0.748, 0.756, 0.756, 0.748, 0.748, 0.760, 0.760, 0.748, 0.760, 0.752, 0.760, 0.752, 0.752, 0.760, 0.752, 0.760, 0.752, 0.760, 0.756, 0.748, 0.748, 0.756, 0.744, 0.752, 0.748, 0.756, 0.740, 0.752, 0.752, 0.740, 0.740, 0.748, 0.744, 0.736, 0.740, 0.732, 0.732, 0.740, 0.736, 0.724, 0.724, 0.732, 0.732, 0.716, 0.716, 0.724, 0.720, 0.708, 0.708, 0.716, 0.716, 0.700, 0.700, 0.708, 0.708, 0.688, 0.692, 0.696, 0.696, 0.680, 0.680, 0.688, 0.684, 0.668, 0.676, 0.668, 0.676, 0.656, 0.656, 0.664, 0.664, 0.644, 0.644, 0.652, 0.648, 0.632, 0.636, 0.628, 0.636, 0.620, 0.624, 0.612, 0.620, 0.612, 0.608, 0.600, 0.608, 0.600, 0.600, 0.584, 0.592, 0.584, 0.588, 0.568, 0.572, 0.576, 0.572, 0.552, 0.564, 0.552, 0.560, 0.540, 0.536, 0.548, 0.540, 0.524, 0.524, 0.532, 0.532, 0.508, 0.508, 0.516, 0.516, 0.492, 0.488, 0.496, 0.500, 0.472, 0.484, 0.472, 0.480, 0.460, 0.464, 0.460, 0.464, 0.448, 0.452, 0.440, 0.448, 0.428, 0.432, 0.424, 0.424, 0.432, 0.416, 0.412, 0.420, 0.408, 0.404, 0.396, 0.400, 0.392, 0.392, 0.380, 0.380, 0.392, 0.380, 0.368, 0.376, 0.364, 0.372, 0.352, 0.356, 0.352, 0.356, 0.340, 0.348, 0.340, 0.344, 0.328, 0.328, 0.336, 0.328, 0.316, 0.312, 0.324, 0.324, 0.304, 0.312, 0.304, 0.312, 0.296, 0.292, 0.300, 0.296, 0.284, 0.288, 0.280, 0.288, 0.276, 0.272, 0.284, 0.280, 0.268, 0.264, 0.272, 0.264, 0.272, 0.268, 0.260, 0.264, 0.256, 0.252, 0.256, 0.260, 0.248, 0.252, 0.244, 0.252, 0.248, 0.252, 0.240, 0.244, 0.248, 0.248, 0.240, 0.240, 0.244, 0.244, 0.236, 0.244, 0.236, 0.236, 0.244, 0.236, 0.240, 0.240, 0.236, 0.232, 0.244, 0.244, 0.236, 0.248, 0.236, 0.236, 0.244, 0.244, 0.240, 0.236, 0.248, 0.248, 0.240, 0.240, 0.248, 0.244, 0.252, 0.252, 0.244, 0.248, 0.260, 0.252, 0.256, 0.256, 0.260, 0.256, 0.264, 0.268, 0.260, 0.260, 0.272, 0.268, 0.276, 0.268, 0.276, 0.288, 0.276, 0.280, 0.284, 0.284, 0.296, 0.288, 0.296, 0.292, 0.304, 0.300, 0.304, 0.300, 0.316, 0.308, 0.316, 0.312, 0.328, 0.320, 0.328, 0.320, 0.340, 0.332, 0.340, 0.332, 0.352, 0.348, 0.344, 0.348, 0.364, 0.360, 0.368, 0.360, 0.376, 0.376, 0.372, 0.376, 0.392, 0.388, 0.396, 0.388, 0.404, 0.412, 0.400, 0.404, 0.420, 0.424, 0.416, 0.416, 0.428, 0.432, 0.440, 0.432, 0.448, 0.448, 0.456, 0.456, 0.448, 0.460, 0.472, 0.472, 0.464, 0.476, 0.488, 0.480, 0.492, 0.488, 0.504, 0.504, 0.496, 0.500, 0.520, 0.512, 0.520, 0.520, 0.532, 0.536, 0.528, 0.532, 0.552, 0.544, 0.552, 0.548, 0.568, 0.560, 0.572, 0.564, 0.584, 0.584, 0.576, 0.576, 0.596, 0.588, 0.600, 0.592, 0.616, 0.616, 0.604, 0.608, 0.624, 0.620, 0.628, 0.620, 0.644, 0.636, 0.640, 0.636, 0.656, 0.656, 0.648, 0.648, 0.664, 0.660, 0.668, 0.660, 0.676, 0.680, 0.672, 0.680, 0.676, 0.680, 0.696, 0.684, 0.692, 0.692, 0.700, 0.696, 0.700, 0.696, 0.712, 0.704, 0.712, 0.704, 0.720, 0.712, 0.724, 0.716, 0.732, 0.724, 0.732, 0.724, 0.740, 0.728, 0.736, 0.728, 0.740, 0.736, 0.744, 0.736, 0.748, 0.748, 0.740, 0.736, 0.752, 0.752, 0.740, 0.752, 0.744, 0.748, 0.756, 0.748, 0.760, 0.760, 0.748, 0.760, 0.748, 0.760, 0.752, 0.760, 0.752, 0.752, 0.756, 0.760, 0.752, 0.760, 0.752, 0.760, 0.752, 0.756, 0.748, 0.756, 0.748, 0.744, 0.756, 0.748, 0.752, 0.752, 0.740, 0.740, 0.752, 0.748, 0.736, 0.748, 0.736, 0.744, 0.732, 0.728, 0.736, 0.740, 0.724, 0.736, 0.724, 0.728, 0.716, 0.724, 0.716, 0.728, 0.712, 0.708, 0.720, 0.716, 0.700, 0.708, 0.700, 0.708, 0.692, 0.688, 0.696, 0.696, 0.676, 0.680, 0.688, 0.684, 0.672, 0.676, 0.664, 0.676, 0.656, 0.656, 0.660, 0.664, 0.652, 0.652, 0.640, 0.640, 0.652, 0.640, 0.628, 0.636, 0.628, 0.628, 0.616, 0.616, 0.624, 0.616, 0.600, 0.600, 0.612, 0.604, 0.584, 0.592, 0.584, 0.596, 0.568, 0.572, 0.580, 0.576, 0.556, 0.564, 0.556, 0.560, 0.540, 0.536, 0.548, 0.544, 0.524, 0.532, 0.520, 0.528, 0.508, 0.516, 0.504, 0.516, 0.492, 0.492, 0.496, 0.500, 0.476, 0.484, 0.476, 0.476, 0.464, 0.460, 0.468, 0.468, 0.452, 0.452, 0.440, 0.452, 0.436, 0.436, 0.428, 0.424, 0.436, 0.424, 0.408, 0.420, 0.408, 0.412, 0.396, 0.404, 0.396, 0.404, 0.384, 0.376, 0.388, 0.384, 0.364, 0.372, 0.368, 0.376, 0.356, 0.352, 0.360, 0.356, 0.340, 0.336, 0.348, 0.344, 0.328, 0.336, 0.324, 0.332, 0.316, 0.324, 0.312, 0.324, 0.304, 0.312, 0.304, 0.308, 0.296, 0.296, 0.292, 0.300, 0.284, 0.292, 0.280, 0.288, 0.280, 0.280, 0.272, 0.280, 0.268, 0.264, 0.276, 0.272, 0.268, 0.268, 0.260, 0.256, 0.264, 0.260, 0.252, 0.260, 0.248, 0.260, 0.248, 0.256, 0.248, 0.252, 0.240, 0.244, 0.248, 0.248, 0.240, 0.248, 0.240, 0.240, 0.248, 0.236, 0.244, 0.240, 0.240, 0.236, 0.244, 0.236, 0.244, 0.240, 0.236, 0.236, 0.244, 0.236, 0.244, 0.236, 0.248, 0.244, 0.236, 0.240, 0.248, 0.244, 0.252, 0.240, 0.248, 0.244, 0.252, 0.244, 0.252, 0.248, 0.256, 0.248, 0.260, 0.252, 0.264, 0.256, 0.264, 0.264, 0.260, 0.260, 0.272, 0.264, 0.272, 0.276, 0.268, 0.280, 0.276, 0.276, 0.284, 0.280, 0.292, 0.284, 0.296, 0.288, 0.304, 0.304, 0.296, 0.300, 0.316, 0.316, 0.308, 0.312, 0.328, 0.324, 0.316, 0.320, 0.336, 0.336, 0.328, 0.332, 0.352, 0.352, 0.344, 0.348, 0.360, 0.364, 0.356, 0.360, 0.376, 0.372, 0.380, 0.372, 0.388, 0.392, 0.384, 0.384, 0.400, 0.400, 0.412, 0.404, 0.412, 0.416, 0.424, 0.416, 0.424, 0.428, 0.436, 0.432, 0.440, 0.440, 0.456, 0.448, 0.456, 0.452, 0.472, 0.464, 0.476, 0.468, 0.488, 0.488, 0.476, 0.480, 0.504, 0.496, 0.504, 0.496, 0.520, 0.520, 0.512, 0.516, 0.536, 0.536, 0.528, 0.532, 0.556, 0.552, 0.544, 0.544, 0.568, 0.560, 0.568, 0.560, 0.584, 0.576, 0.584, 0.576, 0.600, 0.600, 0.592, 0.592, 0.612, 0.608, 0.616, 0.604, 0.624, 0.628, 0.620, 0.620, 0.636, 0.632, 0.644, 0.632, 0.644, 0.644, 0.656, 0.648, 0.656, 0.660, 0.668, 0.668, 0.660, 0.664, 0.680, 0.680, 0.668, 0.676, 0.688, 0.688, 0.684, 0.688, 0.704, 0.696, 0.704, 0.692, 0.712, 0.712, 0.704, 0.704, 0.720, 0.720, 0.712, 0.712, 0.724, 0.720, 0.732, 0.720, 0.736, 0.736, 0.728, 0.728, 0.740, 0.732, 0.744, 0.732, 0.748, 0.740, 0.748, 0.740, 0.752, 0.744, 0.752, 0.744, 0.756, 0.756, 0.744, 0.756, 0.748, 0.760, 0.748, 0.752, 0.756, 0.748, 0.760, 0.752, 0.760, 0.752, 0.760, 0.760, 0.752, 0.752, 0.760, 0.760, 0.748, 0.756, 0.748, 0.748, 0.756, 0.756, 0.744, 0.748, 0.756, 0.756, 0.744, 0.740, 0.752, 0.752, 0.736, 0.740, 0.748, 0.744, 0.732, 0.740, 0.728, 0.736, 0.724, 0.724, 0.736, 0.736, 0.716, 0.724, 0.716, 0.728, 0.712, 0.708, 0.716, 0.716, 0.700, 0.708, 0.700, 0.704, 0.692, 0.688, 0.700, 0.688, 0.696, 0.688, 0.680, 0.688, 0.676, 0.672, 0.664, 0.668, 0.676, 0.656, 0.664, 0.664, 0.656, 0.656, 0.640, 0.640, 0.652, 0.648, 0.628, 0.636, 0.628, 0.632, 0.612, 0.624, 0.616, 0.616, 0.600, 0.608, 0.600, 0.604, 0.588, 0.584, 0.596, 0.588, 0.568, 0.576, 0.568, 0.576, 0.556, 0.564, 0.552, 0.564, 0.540, 0.548, 0.540, 0.548, 0.528, 0.532, 0.524, 0.532, 0.512, 0.516, 0.504, 0.516, 0.500, 0.500, 0.488, 0.496, 0.484, 0.484, 0.472, 0.484, 0.472, 0.468, 0.456, 0.460, 0.468, 0.452, 0.444, 0.440, 0.448, 0.448, 0.424, 0.428, 0.436 };
float w[len];

float* hanning(int N, short itype)
{
        int half, i, idx, n;
        float* w;

        w = (float*)calloc(N, sizeof(float));
        memset(w, 0, N * sizeof(float));

        if (itype == 1)    //periodic function
                n = N - 1;
        else
                n = N;

        if (n % 2 == 0)
        {
                half = n / 2;
                for (i = 0; i < half; i++) //CALC_HANNING   Calculates Hanning window samples.
                        w[i] = 0.5 * (1 - cos(2 * PI * (i + 1) / (n + 1)));

                idx = half - 1;
                for (i = half; i < n; i++) {
                        w[i] = w[idx];
                        idx--;
                }
        }
        else
        {
                half = (n + 1) / 2;
                for (i = 0; i < half; i++) //CALC_HANNING   Calculates Hanning window samples.
                        w[i] = 0.5 * (1 - cos(2 * PI * (i + 1) / (n + 1)));

                idx = half - 2;
                for (i = half; i < n; i++) {
                        w[i] = w[idx];
                        idx--;
                }
        }

        if (itype == 1)    //periodic function
        {
                for (i = N - 1; i >= 1; i--)
                        w[i] = w[i - 1];
                w[0] = 0.0;
        }
        return(w);
}

// Function to calculate the DFT
void calculateDFT()
{
        float Xr[len];
        float Xi[len];
        int i, k, n, N = 1200;

        float x[len];
        double f[len];
        float mean, total_sum = 0, fmax;
        double fdelta;

        i = len;
        while (i-- > 0) {
                total_sum += xn[i];
        }
        mean = total_sum / len;

        i = len;
        while (i-- > 0) {
                x[i] = (xn[i] - mean);
        }

        fmax = (double)(len - 1) / len / 0.5E-6;
        fdelta = fmax / (len - 1);
        i = len;
        f[0] = 0;
        for (i = 1; i < len; i++) {
                f[i] = f[i - 1] + fdelta;
        }

        for (k = 0; k < len; k++) {
                Xr[k] = 0;
                Xi[k] = 0;
                for (n = 0; n < len; n++) {
                        Xr[k]
                                = (Xr[k]
                                        + x[n] * cos(2 * 3.141592 * k * n / N));
                        Xi[k]
                                = (Xi[k]
                                        - x[n] * sin(2 * 3.141592 * k * n / N));
                }

                //printf("(%f) + j(%f)\n", Xr[k], Xi[k]);

                printf("%f %f\n", sqrt(Xr[k] * Xr[k] + Xi[k] * Xi[k]), f[k]);
        }
}


int main()
{
        calculateDFT();

        return 0;
}
 
Nahoru Odpovědět
31.1.2023 2:14
Avatar
DarkCoder
Člen
Avatar
Odpovídá na Caster
DarkCoder:31.1.2023 11:56

Zde máš optimalizovanou symetrickou funkci hannWindow().

#include <math.h>

#define PI 3.14159265358979323846f

void hannWindow(float *window, int size) {
    float argument = 2 * PI / (size - 1);
    for (int i = 0; i < size; i++) {
        window[i] = 0.5f * (1 - cosf(argument * i));
    }
}

A volání funkce:

hannWindow(xn, len);
Nahoru Odpovědět
31.1.2023 11:56
"I ta nejlepší poučka postrádá na významu, není-li patřičně předána." - DarkCoder
Avatar
DarkCoder
Člen
Avatar
Odpovídá na Caster
DarkCoder:31.1.2023 12:04

Asymetrická verze funkce hannWindow() je pak zde:

#include <math.h>

#define PI 3.14159265358979323846f

void hannWindow(float *window, int size, int offset) {
    float argument = 2 * PI / (size - 1);
    for (int i = 0; i < size; i++) {
        window[i] = 0.5f * (1 - cosf(argument * (i + offset)));
    }
}
Nahoru Odpovědět
31.1.2023 12:04
"I ta nejlepší poučka postrádá na významu, není-li patřičně předána." - DarkCoder
Avatar
Caster
Člen
Avatar
Odpovídá na DarkCoder
Caster:31.1.2023 19:57

Díky za pomoc, vyzkouším a dám vědět.

 
Nahoru Odpovědět
31.1.2023 19:57
Avatar
Caster
Člen
Avatar
Caster:4.2.2023 16:11

Ještě mám problém, jak dotatečně pojmenovat strukturu definivanou v dsp.h ve svém programu, nechci přepisovat dsp.h:

typedef struct
{
        int16_t re;           // real portion (a)
        int16_t im;           // imaginary portion (b)
} int16c;                   // Q15 complex number (a + bi)

V programu potřebuji zadat ADC hodnoty signálu jako pole komplexních čísel, kde im = 0. Pro ukázku jsem použil jen prvních 8 hodnot z 1 024.

int16c DataIn[] = {{0x6d9, 0}, {0x7ad, 0}, {0x856, 0}, {0x93c, 0}, {0xa28, 0}, {0xb07, 0}, {0xbd7, 0}, {0xc68, 0}};

int16c *ifftDin;

int main(void) {
    ifftDin = &DataIn;
    DSP_TransformFFT16(ifftDout, ifftDin, ifftc, iscratch, ilog2N);
    ...
}

Při kompilaci ale na řádku s "ifftDin = &DataIn" dostanu chybové hlášení:

../src/main.c:65:13: error: assignment to 'int16c *' {aka 'struct <anonymous> *'} from incompatible pointer type 'int16c (*)[1024]' {aka 'struct <anonymous> (*)[1024]'} [-Werror=incompatible-pointer-types]

Jak kód upravit?

 
Nahoru Odpovědět
4.2.2023 16:11
Avatar
DarkCoder
Člen
Avatar
Odpovídá na Caster
DarkCoder:4.2.2023 17:51

Máš definovaný nový typ, nepotřebuješ znát tag struktury. Chyba říká něco jiného než si myslíš. Chyba není v anonymní struktuře, ale v tom, že dáváš před identifikátor DataIn &, který tam pochopitelně být nemá neboť se jedná o pole.

Nahoru Odpovědět
4.2.2023 17:51
"I ta nejlepší poučka postrádá na významu, není-li patřičně předána." - DarkCoder
Děláme co je v našich silách, aby byly zdejší diskuze co nejkvalitnější. Proto do nich také mohou přispívat pouze registrovaní členové. Pro zapojení do diskuze se přihlas. Pokud ještě nemáš účet, zaregistruj se, je to zdarma.

Zobrazeno 11 zpráv z 11.