PentHertz
/
srsLTE
mirror of https://github.com/PentHertz/srsLTE.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #20 from suttonpd/master 

Adding arbitrary resampler with tests and benchmark
This commit is contained in:
Ismael Gómez-Miguelez 2014-04-17 10:06:02 +01:00
parent 463f0018c5 ae50c990ec
commit fc75cdcddc
9 changed files with 415 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

49
lte/include/lte/resampling/resample_arb.h Normal file
View File

@ -0,0 +1,49 @@
/**
*
* \section COPYRIGHT
*
* Copyright 2013-2014 The libLTE Developers. See the
* COPYRIGHT file at the top-level directory of this distribution.
*
* \section LICENSE
*
* This file is part of the libLTE library.
*
* libLTE is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as
* published by the Free Software Foundation, either version 3 of
* the License, or (at your option) any later version.
*
* libLTE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* A copy of the GNU Lesser General Public License can be found in
* the LICENSE file in the top-level directory of this distribution
* and at http://www.gnu.org/licenses/.
*
*/
#ifndef RESAMPLE_ARB_
#define RESAMPLE_ARB_
#include <stdint.h>
#include <complex.h>
typedef _Complex float cf_t;
#define RESAMPLE_ARB_N 32 // Polyphase filter rows
#define RESAMPLE_ARB_M 8 // Polyphase filter columns
typedef struct {
float rate; // Resample rate
float step; // Step increment through filter
float acc; // Index into filter
cf_t reg[RESAMPLE_ARB_M]; // Our window of samples
}resample_arb_t;
void resample_arb_init(resample_arb_t *q, float rate);
int resample_arb_compute(resample_arb_t *q, cf_t *input, cf_t *output, int n_in);
#endif //RESAMPLE_ARB_

114
lte/lib/resampling/src/resample_arb.c Normal file
View File

@ -0,0 +1,114 @@
/**
*
* \section COPYRIGHT
*
* Copyright 2013-2014 The libLTE Developers. See the
* COPYRIGHT file at the top-level directory of this distribution.
*
* \section LICENSE
*
* This file is part of the libLTE library.
*
* libLTE is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as
* published by the Free Software Foundation, either version 3 of
* the License, or (at your option) any later version.
*
* libLTE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* A copy of the GNU Lesser General Public License can be found in
* the LICENSE file in the top-level directory of this distribution
* and at http://www.gnu.org/licenses/.
*
*/
#include <math.h>
#include <string.h>
#include "lte/resampling/resample_arb.h"
#include "lte/utils/debug.h"
float polyfilt[RESAMPLE_ARB_N][RESAMPLE_ARB_M] =
{{0,0.002400347599485495,-0.006922416132556366,0.0179104136912176,0.99453086623794,-0.008521087756729117,0.0008598969867484128,0.0004992625165376107},
{-0.001903604727400391,0.004479591950094871,-0.01525319260830623,0.04647449496926549,0.9910477342662829,-0.03275243420114668,0.008048813755373533,-0.001216900416836847},
{-0.001750442300940216,0.006728826416921727,-0.02407540632178267,0.07708575473589654,0.9841056525667189,-0.05473739187922162,0.01460652754040275,-0.002745266140572769},
{-0.001807302702047332,0.009130494591071001,-0.03332524241797466,0.1096377743266821,0.9737536341125557,-0.07444712408657775,0.0205085154280773,-0.004077596041064956},
{-0.002005048395707184,0.01166717081713966,-0.04292591596219218,0.1440068251440274,0.9600641782991848,-0.09187282739868929,0.02573649090563353,-0.005218582609802016},
{-0.002303606593778858,0.01431549924598744,-0.05279365431190471,0.1800496557331084,0.9431333663677619,-0.107022659158021,0.03028295984887362,-0.006178495199082676},
{-0.002673824939871474,0.01705308556587308,-0.06283262272105428,0.2176066284621515,0.9230793648715676,-0.1199244901666642,0.03414559686566755,-0.006951034565279722},
{-0.003099294850834212,0.01985022208215894,-0.07294100481214681,0.2564996326404178,0.9000412628278351,-0.1306207218667012,0.03733448898243594,-0.007554128492547957},
{-0.003564415127592977,0.02267702927238637,-0.08300660359189582,0.2965368855416621,0.874178327911914,-0.1391710078036285,0.03986356581130199,-0.007992692892036229},
{-0.004060222849803048,0.02549633298710617,-0.09291211362877161,0.3375102453232687,0.8456688916962831,-0.1456479420005528,0.04175429631970297,-0.00827816851803391},
{-0.004574770330592958,0.02827133859849648,-0.1025308044303699,0.379200979038518,0.8147072909915275,-0.1501397816831265,0.04303284689167361,-0.008423951215857236},
{-0.005100453570725489,0.0309593063669097,-0.1117325201826136,0.4213772298764557,0.7815036335229155,-0.1527444636485677,0.0437349357738538,-0.008441129184587313},
{-0.005625340687997995,0.03351877124912608,-0.1203802123857154,0.463798444519773,0.7462810500374923,-0.1535722712206338,0.04389261707032224,-0.008345136869130621},
{-0.006140888413286479,0.03590276234084221,-0.1283350405000867,0.5062161107254219,0.7092744343987509,-0.152741400055636,0.04355110899623683,-0.008147838953503964},
{-0.006634012711933725,0.03806645580467819,-0.1354549138065051,0.5483763739419955,0.6707272107491931,-0.1503798528838644,0.0427502160096194,-0.007865132034489236},
{-0.007094909626785393,0.03996043743797257,-0.1415969539549883,0.5900207719663111,0.6308907308946271,-0.1466186670140106,0.04153829696698895,-0.007508971112586246},
{-0.007508971112586246,0.04153829696698895,-0.1466186670140106,0.6308907308946271,0.5900207719663111,-0.1415969539549883,0.03996043743797257,-0.007094909626785393},
{-0.007865132034489236,0.0427502160096194,-0.1503798528838644,0.6707272107491931,0.5483763739419955,-0.1354549138065051,0.03806645580467819,-0.006634012711933725},
{-0.008147838953503964,0.04355110899623683,-0.152741400055636,0.7092744343987509,0.5062161107254219,-0.1283350405000867,0.03590276234084221,-0.006140888413286479},
{-0.008345136869130621,0.04389261707032224,-0.1535722712206338,0.7462810500374923,0.463798444519773,-0.1203802123857154,0.03351877124912608,-0.005625340687997995},
{-0.008441129184587313,0.0437349357738538,-0.1527444636485677,0.7815036335229155,0.4213772298764557,-0.1117325201826136,0.0309593063669097,-0.005100453570725489},
{-0.008423951215857236,0.04303284689167361,-0.1501397816831265,0.8147072909915275,0.379200979038518,-0.1025308044303699,0.02827133859849648,-0.004574770330592958},
{-0.00827816851803391,0.04175429631970297,-0.1456479420005528,0.8456688916962831,0.3375102453232687,-0.09291211362877161,0.02549633298710617,-0.004060222849803048},
{-0.007992692892036229,0.03986356581130199,-0.1391710078036285,0.874178327911914,0.2965368855416621,-0.08300660359189582,0.02267702927238637,-0.003564415127592977},
{-0.007554128492547957,0.03733448898243594,-0.1306207218667012,0.9000412628278351,0.2564996326404178,-0.07294100481214681,0.01985022208215894,-0.003099294850834212},
{-0.006951034565279722,0.03414559686566755,-0.1199244901666642,0.9230793648715676,0.2176066284621515,-0.06283262272105428,0.01705308556587308,-0.002673824939871474},
{-0.006178495199082676,0.03028295984887362,-0.107022659158021,0.9431333663677619,0.1800496557331084,-0.05279365431190471,0.01431549924598744,-0.002303606593778858},
{-0.005218582609802016,0.02573649090563353,-0.09187282739868929,0.9600641782991848,0.1440068251440274,-0.04292591596219218,0.01166717081713966,-0.002005048395707184},
{-0.004077596041064956,0.0205085154280773,-0.07444712408657775,0.9737536341125557,0.1096377743266821,-0.03332524241797466,0.009130494591071001,-0.001807302702047332},
{-0.002745266140572769,0.01460652754040275,-0.05473739187922162,0.9841056525667189,0.07708575473589654,-0.02407540632178267,0.006728826416921727,-0.001750442300940216},
{-0.001216900416836847,0.008048813755373533,-0.03275243420114668,0.9910477342662829,0.04647449496926549,-0.01525319260830623,0.004479591950094871,-0.001903604727400391},
{0.0004992625165376107,0.0008598969867484128,-0.008521087756729117,0.99453086623794,0.0179104136912176,-0.006922416132556366,0.002400347599485495,0}};
// TODO: use lte/utils/vector.h and Volk
cf_t dot_prod(cf_t* x, float *y, int len)
{
cf_t res = 0+0*I;
for(int i=0;i<len;i++){
res += x[i]*y[i];
}
return res;
}
// Right-shift our window of samples
void push(resample_arb_t *q, cf_t x)
{
memmove(&q->reg[1], &q->reg[0], (RESAMPLE_ARB_M-1)*sizeof(cf_t));
q->reg[0] = x;
}
// Initialize our struct
void resample_arb_init(resample_arb_t *q, float rate){
memset(q->reg, 0, RESAMPLE_ARB_M*sizeof(cf_t));
q->acc = 0.0;
q->rate = rate;
q->step = (1/rate)*RESAMPLE_ARB_N;
}
// Resample a block of input data
int resample_arb_compute(resample_arb_t *q, cf_t *input, cf_t *output, int n_in){
int cnt = 0;
int n_out = 0;
int idx = 0;
while(cnt < n_in)
{
*output = dot_prod(q->reg, polyfilt[idx], RESAMPLE_ARB_M);
output++;
n_out++;
q->acc += q->step;
idx = (int)roundf(q->acc);
while(idx >= RESAMPLE_ARB_N){
q->acc -= RESAMPLE_ARB_N;
idx -= RESAMPLE_ARB_N;
if(cnt < n_in)
push(q, input[cnt++]);
}
}
return n_out;
}

35
lte/lib/resampling/test/CMakeLists.txt Normal file
View File

@ -0,0 +1,35 @@
#
# Copyright 2012-2013 The libLTE Developers. See the
# COPYRIGHT file at the top-level directory of this distribution.
#
# This file is part of the libLTE library.
#
# libLTE is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# libLTE is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# A copy of the GNU Lesser General Public License can be found in
# the LICENSE file in the top-level directory of this distribution
# and at http://www.gnu.org/licenses/.
#
########################################################################
# RESAMPLING TEST
########################################################################
ADD_EXECUTABLE(resample_arb_test resample_arb_test.c)
TARGET_LINK_LIBRARIES(resample_arb_test lte)
ADD_EXECUTABLE(resample_arb_bench resample_arb_bench.c)
TARGET_LINK_LIBRARIES(resample_arb_bench lte)
ADD_TEST(resample resample_arb_test)

38
lte/lib/resampling/test/resample_arb_bench.c Normal file
View File

@ -0,0 +1,38 @@
#include <stdio.h>
#include <stdlib.h>
#include <strings.h>
#include <unistd.h>
#include <math.h>
#include <time.h>
#include "lte.h"
#include "lte/resampling/resample_arb.h"
typedef _Complex float cf_t;
int main(int argc, char **argv) {
int N=10000000;
float rate = 24.0/25.0;
cf_t *in = malloc(N*sizeof(cf_t));
cf_t *out = malloc(N*sizeof(cf_t));
for(int i=0;i<N;i++)
in[i] = sin(i*2*M_PI/100);
resample_arb_t r;
resample_arb_init(&r, rate);
clock_t start = clock(), diff;
int n_out = resample_arb_compute(&r, in, out, N);
diff = clock() - start;
int msec = diff * 1000 / CLOCKS_PER_SEC;
float thru = (CLOCKS_PER_SEC/(float)diff)*(N/1e6);
printf("Time taken %d seconds %d milliseconds\n", msec/1000, msec%1000);
printf("Rate = %f MS/sec\n", thru);
free(in);
free(out);
printf("Done\n");
exit(0);
}

63
lte/lib/resampling/test/resample_arb_test.c Normal file
View File

@ -0,0 +1,63 @@
#include <stdio.h>
#include <stdlib.h>
#include <strings.h>
#include <unistd.h>
#include <math.h>
#include <time.h>
#include "lte.h"
#include "lte/resampling/resample_arb.h"
typedef _Complex float cf_t;
int main(int argc, char **argv) {
int N = 100; // Number of sinwave samples
int delay = 5; // Delay of our resampling filter
float down = 25.0; // Downsampling rate
for(int up=1;up<down;up++)
{
float rate = up/down;
printf("Testing rate %f\n", rate);
cf_t *in = malloc(N*sizeof(cf_t));
if(!in) {
perror("malloc");
exit(-1);
}
cf_t *out = malloc(N*sizeof(cf_t));
if(!out) {
perror("malloc");
exit(-1);
}
for(int i=0;i<N;i++)
in[i] = sin(i*2*M_PI/N);
// Resample
resample_arb_t r;
resample_arb_init(&r, rate);
int n_out = resample_arb_compute(&r, in, out, N);
// Check interp values
for(int i=delay+1;i<n_out;i++){
float idx = i/rate;
int pre = floor(idx)-delay;
int post = ceil(idx)-delay;
int round = roundf(idx)-delay;
float diff = fabs(creal(in[pre])-creal(in[post]));
float diff2 = fabs(creal(out[i])-creal(in[round]));
if(diff2 > diff && pre != post){
printf("Interpolation failed at index %f", idx);
exit(-1);
}
}
free(in);
free(out);
}
printf("Ok\n");
exit(0);
}

BIN
matlab/resampler/arb_filter.mat Normal file
View File

Binary file not shown.

46
matlab/resampler/arb_resample_linear.m Normal file
View File

@ -0,0 +1,46 @@
function out = arb_resample_linear(sig, rate)
N = 32;
M = 8;
step = (1/rate)*N;
load('arb_filter.mat');
figure;plot(Num);title('Filter impulse response');
[h,w] = freqz(Num, 1);
figure;plot(20*log10(abs(h)));title('Filter freq response');
% Create polyphase partition
poly = reshape(Num, N, M);
% Filter
sig = [zeros(1,(M/2)-1) sig];
k=0;
acc=0;
index=0;
frac=0;
out = [];
while k < length(sig)-M
sig_reg = fliplr(sig(k+1:k+M));
filt_reg1 = poly(index+1, :);
filt_reg2 = poly(mod(index+1,N)+1, :);
res1 = sig_reg*filt_reg1';
res2 = sig_reg*filt_reg2';
if index+1 == 32
res = res1;
else
res = res1 + (res2-res1)*frac;
end
out = [out res];
acc = acc+step;
index = fix(acc);
while index >= N
acc = acc - N;
index = index - N;
k = k+1;
end
frac = abs(acc-index);
end

37
matlab/resampler/arb_resample_nearest.m Normal file
View File

@ -0,0 +1,37 @@
function out = arb_resample_nearest(sig, rate)
N = 32;
M = 8;
step = (1/rate)*N;
load('arb_filter.mat');
figure;plot(Num);title('Filter impulse response');
[h,w] = freqz(Num, 1);
figure;plot(20*log10(abs(h)));title('Filter frequency response');
% Create polyphase partition
poly = reshape(Num, N, M);
% Filter
sig = [zeros(1,(M/2)-1) sig];
k=0;
acc=0;
index=0;
frac=0;
out = [];
while k < length(sig)-M
sig_reg = fliplr(sig(k+1:k+M));
filt_reg1 = poly(index+1, :);
res = sig_reg*filt_reg1';
out = [out res];
acc = acc+step;
index = round(acc);
while index >= N
acc = acc - N;
index = index - N;
k = k+1;
end
end

33
matlab/resampler/arb_resample_test.m Normal file
View File

@ -0,0 +1,33 @@
close all
clear all
up = 24;
down = 25;
d_rate = up/down;
Fs = 100; % Arbitrary sample rate (used for displays)
Fsin = 1;
% Create a sine wave
t = 0:1/Fs:1-1/Fs;
sig = sin(2*pi*t);
out = arb_resample_nearest(sig,d_rate);
% matlab resample for comparison
out2 = resample(sig, up, down);
figure;hold on;title('Ours and matlabs');
l = min(length(out), length(out2));
stem(out(1:l));
stem(out2(1:l), 'r', 'filled');
diff = out2(1:l)-out(1:l);
figure;plot(diff);title('Difference between ours and matlabs');
figure; hold on;title('Original and resampled - no time scaling');
stem(sig);
stem(out, 'r', 'filled');
% Time align and plot
figure;hold on;title('Original and resampled - with time scaling');
stem((1:75)/Fs,real(sig(1:75)));
stem((1:72)/(Fs*d_rate),real(out(1:72)),'r','filled');
xlabel('Time (sec)');ylabel('Signal value');