improved naive convolution

dev
Oona Räisänen 2016-01-11 10:47:16 +02:00
rodzic 9d0d9b3b02
commit 4761816198
4 zmienionych plików z 128 dodań i 78 usunięć

Wyświetl plik

@ -58,13 +58,16 @@ enum WindowType {
};
enum {
KERNEL_LANCZOS2, KERNEL_LANCZOS3, KERNEL_TENT
KERNEL_LANCZOS2, KERNEL_LANCZOS3, KERNEL_TENT, KERNEL_BOX
};
enum eStreamType {
STREAM_TYPE_FILE, STREAM_TYPE_PA, STREAM_TYPE_STDIN
};
enum {
BORDER_REPEAT, BORDER_ZERO
};
typedef struct {
Point pt;

Wyświetl plik

@ -303,49 +303,82 @@ double sinc (double x) {
return (x == 0.0 ? 1 : sin(M_PI*x) / (M_PI*x));
}
namespace kernel {
Wave Lanczos (size_t kernel_len, size_t a) {
Wave kern(kernel_len);
for (size_t i=0; i<kernel_len; i++) {
double x_kern = (1.0*i/(kernel_len-1) - .5)*2*a;
double x_wind = 2.0*i/(kernel_len-1) - 1;
kern[i] = sinc(x_kern) * sinc(x_wind);
}
return kern;
}
Kernel::Kernel(int type) : m_type(type) {
Wave Tent (size_t kernel_len) {
Wave kern(kernel_len);
for (size_t i=0; i<kernel_len; i++) {
double x = 1.0*i/(kernel_len-1);
kern[i] = 1-2*fabs(x-0.5);
}
return kern;
}
double Kernel::at(double x) const {
double val = 0.0;
if (m_type == KERNEL_LANCZOS2) {
int a = 2;
val = (x >= -a && x <= a) ? sinc(x) * sinc(x/a) : 0;
} else if (m_type == KERNEL_LANCZOS3) {
int a = 3;
val = (x >= -a && x <= a) ? sinc(x) * sinc(x/a) : 0;
} else if (m_type == KERNEL_TENT) {
val = (x >= -1 && x <= 1) ? 1.0-std::fabs(x) : 0.0;
} else if (m_type == KERNEL_BOX) {
val = (x >= -1 && x <= 1) ? 0.5 : 0.0;
}
return val;
}
double Kernel::getHalfWidth() const {
if (m_type == KERNEL_LANCZOS2) {
return 2.0;
}
return 1.0;
}
double complexMag (fftw_complex coeff) {
return sqrt(pow(coeff[0],2) + pow(coeff[1],2));
}
double convolveSingle (const Wave& sig, const Kernel& kernel, double x) {
const int signal_len = sig.size();
const int i_dst = std::round(x);
double result = 0;
if (i_dst < 0 || i_dst >= signal_len)
return result;
for (int i_kern=-3; i_kern<=3; i_kern++) {
int i_src = i_dst + i_kern;
double x_kern = i_src - x;
if (i_src >= 0 && i_src < signal_len) {
result += sig.at(i_src) * kernel.at(x_kern);
}
}
return result;
}
Wave convolve (const Wave& sig, const Wave& kernel, bool wrap_around) {
assert (kernel.size() % 2 == 1);
Wave result(sig.size());
const int signal_len = sig.size();
const int kernel_len = kernel.size();
for (size_t i=0; i<sig.size(); i++) {
Wave result(signal_len);
for (size_t i_kern=0; i_kern<kernel.size(); i_kern++) {
int i_new = i - kernel.size()/2 + i_kern;
for (int i_dst=0; i_dst<signal_len; i_dst++) {
for (int i_kern=0; i_kern<kernel_len; i_kern++) {
int i_src = i_dst - kernel_len/2 + i_kern;
if (wrap_around) {
if (i_new < 0)
i_new += result.size();
result[i_new % result.size()] += sig[i] * kernel[i_kern];
if (i_src < 0) {
i_src = signal_len - (abs(i_src) % signal_len);
} else {
i_src = i_src % signal_len;
}
result.at(i_dst) += sig.at(i_src) * kernel[i_kern];
} else {
if (i_new >= 0 && i_new <= int(result.size()-1))
result[i_new] += sig[i] * kernel[i_kern];
if (i_src >= 0 && i_src < signal_len)
result.at(i_dst) += sig.at(i_src) * kernel[i_kern];
}
}
@ -354,50 +387,56 @@ Wave convolve (const Wave& sig, const Wave& kernel, bool wrap_around) {
return result;
}
Wave upsample (const Wave& orig, size_t factor, int kern_type) {
Wave upsample (const Wave& orig, double factor, int kern_type, int border_treatment) {
Wave kern;
if (kern_type == KERNEL_LANCZOS2) {
kern = kernel::Lanczos(factor*2*2 + 1, 2);
} else if (kern_type == KERNEL_LANCZOS3) {
kern = kernel::Lanczos(factor*3*2 + 1, 3);
} else if (kern_type == KERNEL_TENT) {
kern = kernel::Tent(factor*2 + 1);
const Kernel kern(kern_type);
int orig_size = orig.size();
int result_size = std::ceil(orig_size * factor);
Wave result(result_size);
int halfwidth_samples = std::ceil(kern.getHalfWidth() * factor);
for (int i_src=-5; i_src<orig_size+5; i_src++) {
double x_src = i_src + .5;
double val_src = 0.0;
if (i_src < 0 || i_src >= orig_size) {
if (border_treatment == BORDER_REPEAT) {
val_src = orig.at(i_src < 0 ? 0 : orig_size-1);
}
} else {
val_src = orig.at(i_src);
}
if (val_src == 0.0)
continue;
for (int i_dst=factor*i_src-halfwidth_samples-5; i_dst<=factor*i_src+halfwidth_samples+5; i_dst++) {
if (i_dst >= 0 && i_dst < (int)result_size) {
double x_dst = (i_dst + .5) / factor;
result.at(i_dst) += val_src * kern.at(x_dst - x_src);
}
}
}
size_t orig_size = orig.size();
Wave padded(orig_size * factor);
for (size_t i=0; i<orig_size; i++) {
padded[i * factor] = orig[i];
}
padded.insert(padded.begin(), factor-1, 0);
padded.insert(padded.begin(), orig[0]);
padded.push_back(orig[orig_size-1]);
Wave filtered = convolve(padded, kern);
filtered.erase(filtered.begin(), filtered.begin()+factor/2);
filtered.erase(filtered.end()-factor/2, filtered.end());
return filtered;
return result;
}
Wave deriv (const Wave& wave) {
Wave result;
for (size_t i=1; i<wave.size(); i++)
for (int i=1; i<(int)wave.size(); i++)
result.push_back(wave[i] - wave[i-1]);
return result;
}
std::vector<double> peaks (const Wave& wave, size_t n) {
std::vector<double> peaks (const Wave& wave, int n) {
std::vector<std::pair<double,double> > peaks;
for (size_t i=0; i<wave.size(); i++) {
for (int i=0; i<(int)wave.size(); i++) {
double y1 = (i==0 ? wave[0] : wave[i-1]);
double y2 = wave[i];
double y3 = (i==wave.size()-1 ? wave[wave.size()-1] : wave[i+1]);
double y3 = (i==(int)wave.size()-1 ? wave[wave.size()-1] : wave[i+1]);
if ( fabs(y2) >= fabs(y1) && fabs(y2) >= fabs(y3) )
peaks.push_back({ i + gaussianPeak(y1, y2, y3), wave[i]});
peaks.push_back({ gaussianPeak(wave, i, true), wave[i]});
}
std::sort(peaks.begin(), peaks.end(),
[](std::pair<double,double> a, std::pair<double,double> b) {

Wyświetl plik

@ -24,6 +24,16 @@ namespace window {
Wave Gauss (int);
}
class Kernel {
public:
Kernel(int);
double at(double) const;
double getHalfWidth() const;
private:
int m_type;
Wave m_precalc;
};
class DSP {
public:
@ -62,12 +72,13 @@ class DSP {
};
Wave convolve (const Wave&, const Wave&, bool wrap_around=false);
double convolveSingle (const Wave&, const Kernel&, double);
Wave deriv (const Wave&);
Wave peaks (const Wave&, int);
Wave derivPeaks (const Wave&, int);
Wave rms (const Wave&, int);
Wave upsample (const Wave& orig, size_t factor, int kern_type);
double gaussianPeak (double y1, double y2, double y3);
Wave upsample (const Wave& orig, double factor, int kern_type, int border_treatment=BORDER_ZERO);
double gaussianPeak (const Wave& signal, int idx_peak, bool wrap_around=false);
double power (fftw_complex coeff);
double complexMag (fftw_complex coeff);

Wyświetl plik

@ -72,13 +72,13 @@ void decodeYCbCr (uint32_t src, uint8_t *dst) {
double g = (298.082*Y/256) - (100.291*Cb/256) - (208.120*Cr/256) + 135.576;
double b = (298.082*Y/256) + (516.412*Cb/256) - 276.836;
dst[0] = clip(r);
dst[1] = clip(g);
dst[2] = clip(b);
dst[0] = clipToByte(r);
dst[1] = clipToByte(g);
dst[2] = clipToByte(b);
}
void decodeMono (uint32_t src, uint8_t *dst) {
dst[0] = dst[1] = dst[2] = clip((getChannel(src,0) - 15.5) * 1.172); // mmsstv test images
dst[0] = dst[1] = dst[2] = clipToByte((getChannel(src,0) - 15.5) * 1.172); // mmsstv test images
}
@ -91,24 +91,21 @@ Glib::RefPtr<Gdk::Pixbuf> Picture::renderPixbuf(int width) {
img.at(x) = std::vector<uint32_t>(m.num_lines);
}
Wave signal_up =
(upsample_factor > 1 ? upsample(m_video_signal, upsample_factor, KERNEL_TENT) : m_video_signal);
const Kernel sample_kern(KERNEL_LANCZOS2);
for (PixelSample px : m_pixel_grid) {
double signal_t = (px.t/m_drift + m_starts_at) / m_video_dt * upsample_factor;
double val;
if (signal_t < 0 || signal_t >= signal_up.size()-1) {
val = 0;
} else {
val = signal_up[signal_t];
}
double signal_t = (px.t/m_drift + m_starts_at ) / m_video_dt;
double val = convolveSingle(m_video_signal, sample_kern, signal_t);
int x = px.pt.x;
int y = px.pt.y;
int ch = px.ch;
setChannel(img[x][y], ch, clip(round(val*255)));
if (m_progress >= 1.0*y/(m.num_lines-1))
m_has_line.at(y) = true;
setChannel(img[x][y], ch, clipToByte(round(val*255)));
}
/* chroma reconstruction from 4:2:0 */
@ -121,11 +118,11 @@ Glib::RefPtr<Gdk::Pixbuf> Picture::renderPixbuf(int width) {
column_u.push_back(getChannel(img[x][y], 1));
column_v.push_back(getChannel(img[x][y], 2));
}
column_u_filtered = upsample(column_u, 2, KERNEL_TENT);
column_v_filtered = upsample(column_v, 2, KERNEL_TENT);
for (size_t y=0; y < m.num_lines; y++) {
setChannel(img[x][y], 1, clip(column_u_filtered[y+1]));
setChannel(img[x][y], 2, clip(column_v_filtered[y+1]));
column_u_filtered = upsample(column_u, 2, KERNEL_BOX, BORDER_REPEAT);
column_v_filtered = upsample(column_v, 2, KERNEL_BOX, BORDER_REPEAT);
for (int y=0; y < m.num_lines; y++) {
setChannel(img[x][y], 1, clipToByte(column_u_filtered[y]));
setChannel(img[x][y], 2, clipToByte(column_v_filtered[y]));
}
}
}