SDRAngel  4.11.5
Developer docs for <a href="https://github.com/f4exb/sdrangel">SDRangel<\a>, an Open Source Qt5 / OpenGL 3.0+ SDR and signal analyzer frontend to various hardware.
inthalfbandfiltereof.h
Go to the documentation of this file.
1 // Copyright (C) 2018 F4EXB //
3 // written by Edouard Griffiths //
4 // //
5 // Integer half-band FIR based interpolator and decimator //
6 // This is the even/odd double buffer variant. Really useful only when SIMD is //
7 // used //
8 // //
9 // This program is free software; you can redistribute it and/or modify //
10 // it under the terms of the GNU General Public License as published by //
11 // the Free Software Foundation as version 3 of the License, or //
12 // (at your option) any later version. //
13 // //
14 // This program is distributed in the hope that it will be useful, //
15 // but WITHOUT ANY WARRANTY; without even the implied warranty of //
16 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
17 // GNU General Public License V3 for more details. //
18 // //
19 // You should have received a copy of the GNU General Public License //
20 // along with this program. If not, see <http://www.gnu.org/licenses/>. //
22 
23 #ifndef SDRBASE_DSP_INTHALFBANDFILTEREOF_H_
24 #define SDRBASE_DSP_INTHALFBANDFILTEREOF_H_
25 
26 #include <stdint.h>
27 #include <cstdlib>
28 #include "dsp/dsptypes.h"
29 #include "dsp/hbfiltertraits.h"
30 #include "export.h"
31 
32 template<uint32_t HBFilterOrder>
34 public:
36 
37  bool workDecimateCenter(float *x, float *y)
38  {
39  // insert sample into ring-buffer
40  storeSample(*x, *y);
41 
42  switch(m_state)
43  {
44  case 0:
45  // advance write-pointer
47  // next state
48  m_state = 1;
49  // tell caller we don't have a new sample
50  return false;
51 
52  default:
53  // save result
54  doFIR(x, y);
55  // advance write-pointer
57  // next state
58  m_state = 0;
59  // tell caller we have a new sample
60  return true;
61  }
62  }
63 
64  void myDecimate(float x1, float y1, float *x2, float *y2)
65  {
66  storeSample(x1, y1);
68 
69  storeSample(*x2, *y2);
70  doFIR(x2, y2);
72  }
73 
75  void myInterpolateZeroStuffing(float *x1, float *y1, float *x2, float *y2)
76  {
77  storeSample(*x1, *y1);
78  doFIR(x1, y1);
80 
81  storeSample(0, 0);
82  doFIR(x2, y2);
84  }
85 
87  void myInterpolate(float *x1, float *y1, float *x2, float *y2)
88  {
89  // insert sample into ring double buffer
90  m_samples[m_ptr][0] = *x1;
91  m_samples[m_ptr][1] = *y1;
94 
95  // advance pointer
97  m_ptr++;
98  } else {
99  m_ptr = 0;
100  }
101 
102  // first output sample calculated with the middle peak
105 
106  // second sample calculated with the filter
107  doInterpolateFIR(x2, y2);
108  }
109 
110  void myInterpolateInf(float *x1, float *y1, float *x2, float *y2, float *x3, float *y3, float *x4, float *y4)
111  {
112  myInterpolate(x1, y1, x2, y2);
113  myInterpolate(x3, y3, x4, y4);
114  // rotation
115  qint32 x;
116  x = *x1;
117  *x1 = *y1;
118  *y1 = -x;
119  *x2 = -*x2;
120  *y2 = -*y2;
121  x = *x3;
122  *x3 = -*y3;
123  *y3 = x;
124  }
125 
126  void myInterpolateSup(float *x1, float *y1, float *x2, float *y2, float *x3, float *y3, float *x4, float *y4)
127  {
128  myInterpolate(x1, y1, x2, y2);
129  myInterpolate(x3, y3, x4, y4);
130  // rotation
131  qint32 x;
132  x = *x1;
133  *x1 = -*y1;
134  *y1 = x;
135  *x2 = -*x2;
136  *y2 = -*y2;
137  x = *x3;
138  *x3 = *y3;
139  *y3 = -x;
140  }
141 
142 protected:
143  float m_even[2][HBFIRFilterTraits<HBFilterOrder>::hbOrder]; // double buffer technique
144  float m_odd[2][HBFIRFilterTraits<HBFilterOrder>::hbOrder]; // double buffer technique
145  float m_samples[HBFIRFilterTraits<HBFilterOrder>::hbOrder][2]; // double buffer technique
146 
147  int m_ptr;
148  int m_size;
149  int m_state;
150 
151  void storeSample(float x, float y)
152  {
153  if ((m_ptr % 2) == 0)
154  {
155  m_even[0][m_ptr/2] = x;
156  m_even[1][m_ptr/2] = y;
157  m_even[0][m_ptr/2 + m_size] = x;
158  m_even[1][m_ptr/2 + m_size] = y;
159  }
160  else
161  {
162  m_odd[0][m_ptr/2] = x;
163  m_odd[1][m_ptr/2] = y;
164  m_odd[0][m_ptr/2 + m_size] = x;
165  m_odd[1][m_ptr/2 + m_size] = y;
166  }
167  }
168 
170  {
171  m_ptr = m_ptr + 1 < 2*m_size ? m_ptr + 1: 0;
172  }
173 
174  void doFIR(float *x, float *y)
175  {
176  float iAcc = 0;
177  float qAcc = 0;
178 
179  int a = m_ptr/2 + m_size; // tip pointer
180  int b = m_ptr/2 + 1; // tail pointer
181 
182  for (int i = 0; i < HBFIRFilterTraits<HBFilterOrder>::hbOrder / 4; i++)
183  {
184  if ((m_ptr % 2) == 0)
185  {
186  iAcc += (m_even[0][a] + m_even[0][b]) * HBFIRFilterTraits<HBFilterOrder>::hbCoeffsF[i];
187  qAcc += (m_even[1][a] + m_even[1][b]) * HBFIRFilterTraits<HBFilterOrder>::hbCoeffsF[i];
188  }
189  else
190  {
191  iAcc += (m_odd[0][a] + m_odd[0][b]) * HBFIRFilterTraits<HBFilterOrder>::hbCoeffsF[i];
192  qAcc += (m_odd[1][a] + m_odd[1][b]) * HBFIRFilterTraits<HBFilterOrder>::hbCoeffsF[i];
193  }
194 
195  a -= 1;
196  b += 1;
197  }
198 
199  if ((m_ptr % 2) == 0)
200  {
201  iAcc += m_odd[0][m_ptr/2 + m_size/2] * 0.5f;
202  qAcc += m_odd[1][m_ptr/2 + m_size/2] * 0.5f;
203  }
204  else
205  {
206  iAcc += m_even[0][m_ptr/2 + m_size/2 + 1] * 0.5f;
207  qAcc += m_even[1][m_ptr/2 + m_size/2 + 1] * 0.5f;
208  }
209 
210  *x = iAcc; // HB_SHIFT incorrect do not loose the gained bit
211  *y = qAcc;
212  }
213 
214  void doInterpolateFIR(float *x, float *y)
215  {
216  qint32 iAcc = 0;
217  qint32 qAcc = 0;
218 
219  qint16 a = m_ptr;
220  qint16 b = m_ptr + (HBFIRFilterTraits<HBFilterOrder>::hbOrder / 2) - 1;
221 
222  // go through samples in buffer
223  for (int i = 0; i < HBFIRFilterTraits<HBFilterOrder>::hbOrder / 4; i++)
224  {
225  iAcc += (m_samples[a][0] + m_samples[b][0]) * HBFIRFilterTraits<HBFilterOrder>::hbCoeffsF[i];
226  qAcc += (m_samples[a][1] + m_samples[b][1]) * HBFIRFilterTraits<HBFilterOrder>::hbCoeffsF[i];
227  a++;
228  b--;
229  }
230 
231  *x = iAcc * SDR_RX_SCALED;
232  *y = qAcc * SDR_RX_SCALED;
233  }
234 };
235 
236 template<uint32_t HBFilterOrder>
238 {
240 
241  for (int i = 0; i < 2*m_size; i++)
242  {
243  m_even[0][i] = 0.0f;
244  m_even[1][i] = 0.0f;
245  m_odd[0][i] = 0.0f;
246  m_odd[1][i] = 0.0f;
247  m_samples[i][0] = 0.0f;
248  m_samples[i][1] = 0.0f;
249  }
250 
251  m_ptr = 0;
252  m_state = 0;
253 }
254 
255 #endif /* SDRBASE_DSP_INTHALFBANDFILTEREOF_H_ */
float m_odd[2][HBFIRFilterTraits< HBFilterOrder >::hbOrder]
float m_samples[HBFIRFilterTraits< HBFilterOrder >::hbOrder][2]
void myInterpolate(float *x1, float *y1, float *x2, float *y2)
bool workDecimateCenter(float *x, float *y)
float m_even[2][HBFIRFilterTraits< HBFilterOrder >::hbOrder]
void storeSample(float x, float y)
void myInterpolateInf(float *x1, float *y1, float *x2, float *y2, float *x3, float *y3, float *x4, float *y4)
void myInterpolateZeroStuffing(float *x1, float *y1, float *x2, float *y2)
int32_t i
Definition: decimators.h:244
void doInterpolateFIR(float *x, float *y)
void doFIR(float *x, float *y)
void myDecimate(float x1, float y1, float *x2, float *y2)
void myInterpolateSup(float *x1, float *y1, float *x2, float *y2, float *x3, float *y3, float *x4, float *y4)
#define SDR_RX_SCALED
Definition: dsptypes.h:34