# root/SigProc/trunk/SigProc/utility.cpp@5251

Revision 5251, 4.8 KB (checked in by jgaeddert, 6 years ago)

• Property svn:eol-style set to native
Line
1/****************************************************************************
2Copyright 2005,2006 Virginia Polytechnic Institute and State University
3
4This file is part of the OSSIE Signal Processing Library.
5
6OSSIE Signal Processing Library is free software; you can redistribute it
7and/or modify it under the terms of the GNU General Public License as
9or (at your option) any later version.
10
11OSSIE Signal Processing Library is distributed in the hope that it will be
12useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
13MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14GNU General Public License for more details.
15
16You should have received a copy of the GNU General Public License along
17with OSSIE Signal Processing Library; if not, write to the Free Software
18Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19
20
21****************************************************************************/
22
23#include "SigProc.h"
24
25namespace SigProc {
26
27//-----------------------------------------------------------------------------
28//
29// Dot product definitions
30//
31//-----------------------------------------------------------------------------
32/** \brief Calculate dot product between two floating point arrays
33
34    \f[ z = x \otimes y \f]
35 */
36void dot_product(float *x, float *y, unsigned int N, float &z)
37{
38    z = 0.0f;
39    for (unsigned int i=0; i<N; i++)
40        z += x[i]*y[i];
41}
42
43/** \brief Calculate dot product between two floating short integer arrays
44
45    \f[ z = x \otimes y \f]
46 */
47void dot_product(short *x, short *y, unsigned int N, short &z)
48{
49    float z_tmp(0.0f);
50    for (unsigned int i=0; i<N; i++)
51        z_tmp += float(x[i])*float(y[i]);
52    if ( z_tmp > SHRT_MAX )
53        z = SHRT_MAX;
54    else
55        z = short(z_tmp);
56
57}
58
59
60/** \brief Calculate dot product between a float and a short array
61
62    \f[ z = x \otimes y \f]
63 */
64void dot_product(float *x, short *y, unsigned int N, short &z)
65{
66    float z_tmp;
67
68    dot_product(x, y, N, z_tmp);
69
70    if ( z_tmp > SHRT_MAX )
71        z = SHRT_MAX;
72    else if ( z_tmp < SHRT_MIN )
73        z = SHRT_MIN;
74    else
75        z = short(z_tmp);
76}
77
78
79/** \brief Calculate dot product between a float and a short array
80
81    \f[ z = x \otimes y \f]
82 */
83void dot_product(float *x, short *y, unsigned int N, float &z)
84{
85    z = 0.0f;
86    for (unsigned int i=0; i<N; i++)
87        z += x[i]*float(y[i]);
88}
89
90
91
92//-----------------------------------------------------------------------------
93//
94// Trigonometric functions
95//
96//-----------------------------------------------------------------------------
97/** \brief Calculate inverse tangent
98
99    \f[ \theta = \arctan(y/x);\ \ 0\le\theta< 2\pi \f]
100 */
101void arctan(float &x, float &y, float &theta)
102{
103    if (x==0.0f) {
104        theta = ( y < 0 ) ? 3*PI/2 : PI/2;
105        return;
106    }
107
108    theta = atan(y/x);
109
110    // rotate angle to positive value
111    if (x>0 && y>=0) {
112        // Q1: do nothing
113    } else if (x<0 && y>=0) {
114        // Q2
115        theta += PI;
116    } else if (x<0 && y<=0) {
117        // Q3
118        theta += PI;
119    } else {
120        // Q4
121        theta += TWO_PI;
122    }
123}
124
125
126void rotate(short I_in, short Q_in, float theta, short *I_out, short *Q_out)
127{
128    ///\todo: use CORDIC mixer to perform this operation more efficiently
129    float c = cosf(theta);
130    float s = sinf(theta);
131
132    *I_out = (short)  ( (float) (I_in*c) - (float) (Q_in*s)  );
133    *Q_out = (short)  ( (float) (I_in*s) + (float) (Q_in*c)  );
134}
135
136
137//-----------------------------------------------------------------------------
138//
139// Random number generators
140//
141//-----------------------------------------------------------------------------
142
143float randf() {
144    float x = (float) rand();
145    return x / (float) RAND_MAX;
146}
147
148void randnf(float * i, float * q)
149{
150    // generate two uniform random numbers
151    float u1, u2;
152
153    // ensure u1 does not equal zero
154    do {
155        u1 = randf();
156    } while (u1 == 0.0f);
157
158    u2 = randf();
159
160    float x = sqrt(-2*logf(u1));
161    *i = x * sinf(2*PI*u2);
162    *q = x * cosf(2*PI*u2);
163}
164
165
166
167dump_data::dump_data(const char *filename, long _start_sample, long _number_of_samples) : start_sample(_start_sample), stop_sample(_start_sample + _number_of_samples), current_sample(0)
168{
169
170    out_file = new std::ofstream(filename);
171
172}
173
174dump_data::~dump_data()
175{
176    delete out_file;
177}
178
179
180void dump_data::write_data(float data, const char *msg)
181{
182    if ((current_sample < stop_sample) && (current_sample > start_sample)) {
183        *out_file << msg << data << std::endl;
184    }
185    ++current_sample;
186}
187
188void dump_data::write_data(float a, float b, const char *msg)
189{
190    if ((current_sample < stop_sample) && (current_sample > start_sample)) {
191        *out_file << msg << a << "  " << b << std::endl;
192    }
193    ++current_sample;
194}
195
196} // namespace SigProc
Note: See TracBrowser for help on using the browser.