summaryrefslogtreecommitdiff
path: root/source/application/main/PlatformMath.cc
blob: 00436352f60856cd84497a61e9b07332bac77d49 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
/*
 * Copyright (c) 2021 Arm Limited. All rights reserved.
 * SPDX-License-Identifier: Apache-2.0
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
#include "PlatformMath.hpp"

#if 0 == ARM_DSP_AVAILABLE
    #include <cmath>
    #include <numeric>
#endif /* 0 == ARM_DSP_AVAILABLE */

namespace arm {
namespace app {
namespace math {

    float MathUtils::CosineF32(float radians)
    {
#if ARM_DSP_AVAILABLE
        return arm_cos_f32(radians);
#else /* ARM_DSP_AVAILABLE */
        return cos(radians);
#endif /* ARM_DSP_AVAILABLE */
    }

    float MathUtils::SqrtF32(float input)
    {
#if ARM_DSP_AVAILABLE
        float output = 0.f;
        arm_sqrt_f32(input, &output);
        return output;
#else /* ARM_DSP_AVAILABLE */
        return sqrtf(input);
#endif /* ARM_DSP_AVAILABLE */
    }

    float MathUtils::MeanF32(float* ptrSrc, const uint32_t srcLen)
    {
        if (!srcLen) {
            return 0.f;
        }

#if ARM_DSP_AVAILABLE
        float result = 0.f;
        arm_mean_f32(ptrSrc, srcLen, &result);
        return result;
#else /* ARM_DSP_AVAILABLE */
        float acc = std::accumulate(ptrSrc, ptrSrc + srcLen, 0.0);
        return acc/srcLen;
#endif /* ARM_DSP_AVAILABLE */
    }

    float MathUtils::StdDevF32(float* ptrSrc, const uint32_t srcLen,
                               const float mean)
    {
        if (!srcLen) {
            return 0.f;
        }
#if ARM_DSP_AVAILABLE
        /**
         * Note Standard deviation calculation can be off
         * by > 0.01 but less than < 0.1, according to
         * preliminary findings.
         **/
        UNUSED(mean);
        float stdDev = 0;
        arm_std_f32(ptrSrc, srcLen, &stdDev);
        return stdDev;
#else /* ARM_DSP_AVAILABLE */
        auto VarianceFunction = [=](float acc, const float value) {
            return acc + (((value - mean) * (value - mean))/ srcLen);
        };

        float acc = std::accumulate(ptrSrc, ptrSrc + srcLen, 0.0,
                                    VarianceFunction);

        return sqrtf(acc);
#endif /* ARM_DSP_AVAILABLE */
    }

    bool MathUtils::FftInitF32(const uint16_t fftLen, arm::app::math::FftInstance& fftInstance)
    {
#if ARM_DSP_AVAILABLE
        if (!fftInstance.initialised) {
            arm_status status = arm_rfft_fast_init_f32(&fftInstance.instance, fftLen);

            if (ARM_MATH_SUCCESS != status) {
                return false;
            }
            fftInstance.initialised = true;
        }
#else
        UNUSED(fftLen);
        UNUSED(fftInstance);
#endif /* ARM_DSP_AVAILABLE */
        return true;
    }

    void MathUtils::FftF32(std::vector<float>& input,
                           std::vector<float>& fftOutput,
                           arm::app::math::FftInstance& fftInstance)
    {
#if ARM_DSP_AVAILABLE
        arm_rfft_fast_f32(&fftInstance.instance, input.data(), fftOutput.data(), 0);
#else
        UNUSED(fftInstance);
        const int inputLength = input.size();

        for (int k = 0; k <= inputLength / 2; k++) {
            float sumReal = 0, sumImag = 0;

            for (int t = 0; t < inputLength; t++) {
                auto angle = static_cast<float>(2 * M_PI * t * k / inputLength);
                sumReal += input[t] * cosf(angle);
                sumImag += -input[t] * sinf(angle);
            }

            /* Arrange output to [real0, realN/2, real1, im1, real2, im2, ...] */
            if (k == 0) {
                fftOutput[0] = sumReal;
            } else if (k == inputLength / 2) {
                fftOutput[1] = sumReal;
            } else {
                fftOutput[k*2] = sumReal;
                fftOutput[k*2 + 1] = sumImag;
            };
        }
#endif /* ARM_DSP_AVAILABLE */
    }

    void MathUtils::VecLogarithmF32(std::vector <float>& input,
                                    std::vector <float>& output)
    {
#if ARM_DSP_AVAILABLE
        arm_vlog_f32(input.data(), output.data(),
                     output.size());
#else /* ARM_DSP_AVAILABLE */
        for (auto in = input.begin(), out = output.begin();
             in != input.end() && out != output.end(); ++in, ++out) {
            *out = logf(*in);
        }
#endif /* ARM_DSP_AVAILABLE */
    }

    float MathUtils::DotProductF32(float* srcPtrA, float* srcPtrB,
                                   const uint32_t srcLen)
    {
        float output = 0.f;

#if ARM_DSP_AVAILABLE
        arm_dot_prod_f32(srcPtrA, srcPtrB, srcLen, &output);
#else /* ARM_DSP_AVAILABLE */
        for (uint32_t i = 0; i < srcLen; ++i) {
            output += *srcPtrA++ * *srcPtrB++;
        }
#endif /* ARM_DSP_AVAILABLE */

        return output;
    }

    bool MathUtils::ComplexMagnitudeSquaredF32(float* ptrSrc,
                                               const uint32_t srcLen,
                                               float* ptrDst,
                                               const uint32_t dstLen)
    {
        if (dstLen < srcLen/2) {
            printf_err("dstLen must be greater than srcLen/2");
            return false;
        }

#if ARM_DSP_AVAILABLE
        arm_cmplx_mag_squared_f32(ptrSrc, ptrDst, srcLen/2);
#else /* ARM_DSP_AVAILABLE */
        for (uint32_t j = 0; j < srcLen/2; ++j) {
            const float real = *ptrSrc++;
            const float im = *ptrSrc++;
            *ptrDst++ = real*real + im*im;
        }
#endif /* ARM_DSP_AVAILABLE */
        return true;
    }

} /* namespace math */
} /* namespace app */
} /* namespace arm */