From 6b040e853f28da50027ef403d520bf5d057628f1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ahmet=20=C4=B0nan?= Date: Mon, 5 Sep 2011 23:39:05 +0200 Subject: [PATCH] Initial commit --- Makefile | 25 ++ decode.c | 769 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ encode.c | 355 +++++++++++++++++++++++++ smpte.ppm | Bin 0 -> 230416 bytes 4 files changed, 1149 insertions(+) create mode 100644 Makefile create mode 100644 decode.c create mode 100644 encode.c create mode 100644 smpte.ppm diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..cc398c6 --- /dev/null +++ b/Makefile @@ -0,0 +1,25 @@ + +CC = gcc +CFLAGS = -D_GNU_SOURCE= -W -Wall -O3 -std=c99 -lm -ffast-math + +all: encode decode + ./encode smpte.ppm 8000.wav 8000 + ./encode smpte.ppm 11025.wav 11025 + ./encode smpte.ppm 40000.wav 40000 + ./encode smpte.ppm 44100.wav 44100 + ./encode smpte.ppm 48000.wav 48000 + ./decode 8000.wav 8000.ppm + ./decode 11025.wav 11025.ppm + ./decode 40000.wav 40000.ppm + ./decode 44100.wav 44100.ppm + ./decode 48000.wav 48000.ppm + +clean: + rm -f encode decode {8000,11025,40000,44100,48000}.{ppm,wav} + +encode: encode.c Makefile + $(CC) -o $@ $< $(CFLAGS) + +decode: decode.c Makefile + $(CC) -o $@ $< $(CFLAGS) -DDN=1 -DUP=1 + diff --git a/decode.c b/decode.c new file mode 100644 index 0000000..5485adf --- /dev/null +++ b/decode.c @@ -0,0 +1,769 @@ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +float lerp(float a, float b, float x) +{ + return a - a * x + b * x; +} + +int64_t gcd(int64_t a, int64_t b) +{ + int64_t c; + while ((c = a % b)) { + a = b; + b = c; + } + return b; +} + +float limit(float min, float max, float x) +{ + float tmp = x < min ? min : x; + return tmp > max ? max : tmp; +} + +float sinc(float x) +{ + return 0 == x ? 1.0 : sinf(M_PI * x) / (M_PI * x); +} +float hann(float n, float N) +{ + return 0.5 * (1.0 - cosf(2.0 * M_PI * n / (N - 1.0))); +} +float hamming(float n, float N) +{ + return 0.54 - 0.46 * cosf(2.0 * M_PI * n / (N - 1.0)); +} +float lanczos(float n, float N) +{ + return sinc(2.0 * n / (N - 1.0) - 1.0); +} +float gauss(float n, float N) +{ + float o = 0.35; + return expf(- 1.0/2.0 * powf((n - (N - 1.0) / 2.0) / (o * (N - 1.0) / 2.0), 2.0)); +} + +float i0f(float x) +{ + // converges for -3*M_PI:3*M_PI in less than 20 iterations + float sum = 1.0, val = 1.0, c = 0.0; + for (int n = 1; n < 20; n++) { + float tmp = x / (2 * n); + val *= tmp * tmp; + float y = val - c; + float t = sum + y; + c = (t - sum) - y; + sum = t; + } + return sum; +} +float kaiser(float n, float N) +{ + float a = 2.0; + return i0f(M_PI * a * sqrtf(1.0 - powf((2.0 * n) / (N - 1.0) - 1.0, 2.0))) / i0f(M_PI * a); +} + +typedef struct { + float complex *b; + float *s; + float complex osc; + float complex d; + int offset; + int skip; + int last; + int taps; + int samples; + int L; + int M; +} ddc_t; + +void do_ddc(ddc_t *ddc, float *input, float complex *output) +{ + int in = 0; + ddc->s[ddc->last] = input[in++]; + ddc->last = (ddc->last + 1) < ddc->samples ? ddc->last + 1 : 0; + ddc->skip += ddc->L; + // this works only for L <= M + for (int k = 0; k < ddc->L; k++) { + float complex sum = 0.0; + for (int i = ddc->offset, j = ddc->last; i < ddc->taps; i += ddc->L) { + sum += ddc->b[i] * ddc->s[j]; + j += j ? - 1 : ddc->samples - 1; + } + + ddc->offset = (ddc->offset + ddc->M) % ddc->L; + + while (ddc->skip < ddc->M) { + ddc->s[ddc->last] = input[in++]; + ddc->last = (ddc->last + 1) < ddc->samples ? ddc->last + 1 : 0; + ddc->skip += ddc->L; + } + + ddc->skip %= ddc->M; + output[k] = ddc->osc * sum; + ddc->osc *= ddc->d; +// ddc->osc /= cabsf(ddc->osc); // not really needed + } +} +ddc_t *alloc_ddc(float freq, float bw, float step, int taps, int L, int M, float (*window)(float, float)) +{ + float lstep = step / (float)L; + float ostep = step * (float)M / (float)L; + ddc_t *ddc = malloc(sizeof(ddc_t)); + ddc->taps = taps; + ddc->samples = (taps + L - 1) / L; + ddc->b = malloc(sizeof(float complex) * ddc->taps); + ddc->s = malloc(sizeof(float) * ddc->samples); + ddc->osc = I; + ddc->d = cexpf(-I * 2.0 * M_PI * freq * ostep); + ddc->offset = 0; + ddc->last = 0; + ddc->skip = 0; + ddc->L = L; + ddc->M = M; + for (int i = 0; i < ddc->samples; i++) + ddc->s[i] = 0.0; + float sum = 0.0; + for (int i = 0; i < ddc->taps; i++) { + float N = (float)ddc->taps; + float n = (float)i; + float x = n - (N - 1.0) / 2.0; + float l = 2.0 * M_PI * bw * lstep; + float w = window(n, ddc->taps); + float h = 0.0 == x ? l / M_PI : sinf(l * x) / (x * M_PI); + float b = w * h; + sum += b; + complex float o = cexpf(I * 2.0 * M_PI * freq * lstep * n); + ddc->b[i] = b * o * (float)L; + } + for (int i = 0; i < ddc->taps; i++) + ddc->b[i] /= sum; + return ddc; +} +void free_ddc(ddc_t *ddc) +{ + free(ddc->b); + free(ddc->s); + free(ddc); +} + +typedef struct { + float *s; + int last; + int len; +} delay_t; + +float do_delay(delay_t *d, float input) +{ + d->s[d->last] = input; + d->last = (d->last + 1) < d->len ? d->last + 1 : 0; + return d->s[d->last]; +} + +delay_t *alloc_delay(int samples) +{ + int len = samples + 1; + delay_t *d = malloc(sizeof(delay_t)); + d->s = malloc(sizeof(float) * len); + d->last = 0; + d->len = len; + for (int i = 0; i < len; i++) + d->s[i] = 0.0; + return d; +} +void free_delay(delay_t *delay) +{ + free(delay->s); + free(delay); +} + +void *mmap_file_ro(char *name, size_t *size) +{ + *size = 0; + int fd = open(name, O_RDONLY); + if (fd == -1) { + perror("open"); + return 0; + } + + struct stat sb; + if (fstat(fd, &sb) == -1) { + perror("fstat"); + return 0; + } + + if (!S_ISREG(sb.st_mode)) { + fprintf(stderr, "%s not a file\n", name); + return 0; + } + + void *p = mmap(0, sb.st_size, PROT_READ, MAP_SHARED, fd, 0); + if (p == MAP_FAILED) { + perror("mmap"); + return 0; + } + + if (close(fd) == -1) { + perror ("close"); + return 0; + } + *size = sb.st_size; + return p; +} + +void *mmap_file_rw(char *name, size_t size) +{ + int fd = open(name, O_RDWR|O_CREAT|O_TRUNC, S_IRUSR|S_IWUSR|S_IRGRP|S_IROTH); + if (fd == -1) { + perror("open"); + return 0; + } + + struct stat sb; + if (fstat(fd, &sb) == -1) { + perror("fstat"); + return 0; + } + + if (!S_ISREG(sb.st_mode)) { + fprintf(stderr, "%s not a file\n", name); + return 0; + } + + if (lseek(fd, size - 1, SEEK_SET) == -1) { + perror("lseek"); + return 0; + } + + if (write(fd, "", 1) != 1) { + perror("write"); + return 0; + } + + void *p = mmap(0, size, PROT_READ|PROT_WRITE, MAP_SHARED, fd, 0); + if (p == MAP_FAILED) { + perror("mmap"); + return 0; + } + + if (close(fd) == -1) { + perror ("close"); + return 0; + } + return p; +} +int munmap_file(void *p, size_t size) +{ + if (munmap(p, size) == -1) { + perror("munmap"); + return 0; + } + return 1; +} + +typedef struct { + uint32_t ChunkID; + uint32_t ChunkSize; + uint32_t Format; + uint32_t Subchunk1ID; + uint32_t Subchunk1Size; + uint16_t AudioFormat; + uint16_t NumChannels; + uint32_t SampleRate; + uint32_t ByteRate; + uint16_t BlockAlign; + uint16_t BitsPerSample; + uint32_t Subchunk2ID; + uint32_t Subchunk2Size; +} wav_t; + +uint8_t R_YUV(uint8_t Y, uint8_t U, uint8_t V) +{ + (void)U; + return limit(0.0, 255.0, 0.003906 * ((298.082 * (Y - 16.0)) + (408.583 * (V - 128)))); +} +uint8_t G_YUV(uint8_t Y, uint8_t U, uint8_t V) +{ + return limit(0.0, 255.0, 0.003906 * ((298.082 * (Y - 16.0)) + (-100.291 * (U - 128)) + (-208.12 * (V - 128)))); +} +uint8_t B_YUV(uint8_t Y, uint8_t U, uint8_t V) +{ + (void)V; + return limit(0.0, 255.0, 0.003906 * ((298.082 * (Y - 16.0)) + (516.411 * (U - 128)))); +} + +void process_line(uint8_t *pixel, uint8_t *y_pixel, uint8_t *uv_pixel, int y_width, int uv_width, int width, int height, int n) +{ + // we only process after 2 full lines: on odd lines + if (n % 2) + for (int y = n-1, l = 0; l < 2 && y < height; l++, y++) { + for (int x = 0; x < width; x++) { +#if DN && UP + uint8_t Y = y_pixel[x + l*y_width]; + uint8_t U = uv_pixel[x/2 + uv_width]; + uint8_t V = uv_pixel[x/2]; +#else + float y_xf = (float)x * (float)y_width / (float)width; + float uv_xf = (float)x * (float)uv_width / (float)width; + int y_x0 = y_xf; + int uv_x0 = uv_xf; + int y_x1 = limit(0, y_width, y_xf + 1); + int uv_x1 = limit(0, uv_width, uv_xf + 1); + uint8_t Y = lerp(y_pixel[y_x0 + l*y_width], y_pixel[y_x1 + l*y_width], y_xf - (float)y_x0); + uint8_t U = lerp(uv_pixel[uv_x0 + uv_width], uv_pixel[uv_x1 + uv_width], uv_xf - (float)uv_x0); + uint8_t V = lerp(uv_pixel[uv_x0], uv_pixel[uv_x1], uv_xf - (float)uv_x0); +#endif + uint8_t *p = pixel + 3 * width * y + 3 * x; + p[0] = R_YUV(Y, U, V); + p[1] = G_YUV(Y, U, V); + p[2] = B_YUV(Y, U, V); + } + } +} + +int main(int argc, char **argv) +{ + if (argc != 3) { + fprintf(stderr, "usage: %s \n", argv[0]); + return 1; + } + + size_t wav_size; + char *wav_p = mmap_file_ro(argv[1], &wav_size); + if (!wav_p) + return 1; + + wav_t *wav = (wav_t *)wav_p; + +#if 0 + fprintf(stderr, "\n"); + fprintf(stderr, "ChunkID = 0x%x\n", wav->ChunkID); + fprintf(stderr, "ChunkSize = %d\n", wav->ChunkSize); + fprintf(stderr, "Format = 0x%x\n", wav->Format); + fprintf(stderr, "Subchunk1ID = 0x%x\n", wav->Subchunk1ID); + fprintf(stderr, "Subchunk1Size = %d\n", wav->Subchunk1Size); + fprintf(stderr, "AudioFormat = %d\n", wav->AudioFormat); + fprintf(stderr, "NumChannels = %d\n", wav->NumChannels); + fprintf(stderr, "SampleRate = %d\n", wav->SampleRate); + fprintf(stderr, "ByteRate = %d\n", wav->ByteRate); + fprintf(stderr, "BlockAlign = %d\n", wav->BlockAlign); + fprintf(stderr, "BitsPerSample = %d\n", wav->BitsPerSample); + fprintf(stderr, "Subchunk2ID = 0x%x\n", wav->Subchunk2ID); + fprintf(stderr, "Subchunk2Size = %d\n", wav->Subchunk2Size); + fprintf(stderr, "\n"); +#endif + + if (wav->ChunkID != 0x46464952 || wav->Format != 0x45564157 || + wav->Subchunk1ID != 0x20746d66 || wav->Subchunk1Size != 16 || + wav->AudioFormat != 1 || wav->Subchunk2ID != 0x61746164) { + fprintf(stderr, "unsupported WAV file!\n"); + return 1; + } + if (wav->BitsPerSample != 16) { + fprintf(stderr, "only 16bit WAV supported!\n"); + return 1; + } + + if (wav->NumChannels != 1) { + fprintf(stderr, "only Mono WAV supported!\n"); + return 1; + } + + int samples = wav->Subchunk2Size / 2; + + float rate = wav->SampleRate; + if (rate * 0.088 < 320.0) { + fprintf(stderr, "%.0fhz samplerate too low\n", rate); + return 1; + } + fprintf(stderr, "%.0fhz samplerate\n", rate); + + short *b = (short *)(wav_p + sizeof(wav_t)); + + const float step = 1.0 / rate; + float complex cnt_last = -I; + float complex dat_last = -I; + + float cal_avg = 1900.0; + + int begin_vis_ss = 0; + int begin_vis_lo = 0; + int begin_vis_hi = 0; + int begin_hor_sync = 0; + int begin_cal_break = 0; + int begin_cal_leader = 0; + int begin_sep_evn = 0; + int begin_sep_odd = 0; + int latch_sync = 0; + + const float vis_len = 0.03; + const float hor_sync_len = 0.009; + const float cal_break_len = 0.01; + const float cal_leader_len = 0.3; + const float seperator_len = 0.0045; + int cal_ticks = 0; + int got_cal_break = 0; + int vis_mode = 0; + int dat_mode = 0; + int vis_ticks = 0; + int vis_bit = -1; + int vis_byte = 0; + + int y = 0; + int odd = 0; + int odd_count = 0; + int evn_count = 0; + int first_hor_sync = 0; + +#if DN && UP + // 320 / 0.088 = 160 / 0.044 = 40000 / 11 = 3636.(36)~ pixels per second for Y, U and V + int64_t factor_L = 40000; + int64_t factor_M = 11 * rate; + int64_t factor_D = gcd(factor_L, factor_M); + factor_L /= factor_D; + factor_M /= factor_D; +#endif +#if DN && !UP + int64_t factor_L = 1; + // factor_M * step should be smaller than pixel length + int64_t factor_M = rate * 0.088 / 320.0; +#endif +#if !DN + int64_t factor_L = 1; + int64_t factor_M = 1; +#endif + + // we want odd number of taps, 4 and 2 ms window length gives best results + int cnt_taps = 1 | (int)(rate * factor_L * 0.004); + int dat_taps = 1 | (int)(rate * factor_L * 0.002); + fprintf(stderr, "using %d and %d tap filter\n", cnt_taps, dat_taps); + float drate = rate * (float)factor_L / (float)factor_M; + float dstep = 1.0 / drate; + fprintf(stderr, "using factor of %ld/%ld, working at %.2fhz\n", factor_L, factor_M, drate); + float *cnt_amp = malloc(sizeof(float) * factor_M); + float *dat_amp = malloc(sizeof(float) * factor_M); + float complex *cnt_q = malloc(sizeof(float complex) * factor_L); + float complex *dat_q = malloc(sizeof(float complex) * factor_L); + // same factor to keep life simple and have accurate horizontal sync + ddc_t *cnt_ddc = alloc_ddc(1200.0, 200.0, step, cnt_taps, factor_L, factor_M, kaiser); + ddc_t *dat_ddc = alloc_ddc(1900.0, 800.0, step, dat_taps, factor_L, factor_M, kaiser); + // delay input by phase shift of other filter to synchronize outputs + delay_t *cnt_delay = alloc_delay((dat_taps - 1) / (2 * factor_L)); + delay_t *dat_delay = alloc_delay((cnt_taps - 1) / (2 * factor_L)); + + const float sync_porch_len = 0.003; + const float porch_len = 0.0015; (void)porch_len; + const float y_len = 0.088; + const float uv_len = 0.044; + const float hor_len = 0.15; + int missing_sync = 0; + int seperator_correction = 0; + +#if DEBUG + const int width = (0.150 + 3.0 * sync_porch_len) * drate + 20; + const int height = 256; +#else + const int width = 320; + const int height = 240; +#endif + char ppm_head[32]; + snprintf(ppm_head, 32, "P6 %d %d 255\n", width, height); + size_t ppm_size = strlen(ppm_head) + width * height * 3; + char *ppm_p = mmap_file_rw(argv[2], ppm_size); + memcpy(ppm_p, ppm_head, strlen(ppm_head)); + uint8_t *pixel = (uint8_t *)ppm_p + strlen(ppm_head); + memset(pixel, 0, width * height * 3); + + int hor_ticks = 0; + int y_pixel_x = 0; + int uv_pixel_x = 0; + int y_width = drate * y_len; + int uv_width = drate * uv_len; + uint8_t *y_pixel = malloc(y_width * 2); + memset(y_pixel, 0, y_width * 2); + uint8_t *uv_pixel = malloc(uv_width * 2); + memset(uv_pixel, 0, uv_width * 2); + + for (int ticks = 0, in = 0, out = factor_L; in < (samples + 1 - factor_M); out++, ticks++, hor_ticks++, cal_ticks++, vis_ticks++) { + if (out >= factor_L) { + out = 0; + for (int j = 0; j < factor_M; j++) { + float amp = (float)b[in++] / 32767.0; + cnt_amp[j] = do_delay(cnt_delay, amp); + dat_amp[j] = do_delay(dat_delay, amp); + } + do_ddc(cnt_ddc, cnt_amp, cnt_q); + do_ddc(dat_ddc, dat_amp, dat_q); + } + + float cnt_freq = limit(1100.0, 1300.0, 1200.0 + cargf(cnt_q[out] * conjf(cnt_last)) / (2.0 * M_PI * dstep)); + float dat_freq = limit(1500.0, 2300.0, 1900.0 + cargf(dat_q[out] * conjf(dat_last)) / (2.0 * M_PI * dstep)); + + cnt_last = cnt_q[out]; + dat_last = dat_q[out]; + + const float cal_a = 0.05; + cal_avg = cal_a * dat_freq + (1.0 - cal_a) * cal_avg; + + begin_vis_ss = fabsf(cnt_freq - 1200.0) < 50.0 ? begin_vis_ss + 1 : 0; + begin_vis_lo = fabsf(cnt_freq - 1300.0) < 50.0 ? begin_vis_lo + 1 : 0; + begin_vis_hi = fabsf(cnt_freq - 1100.0) < 50.0 ? begin_vis_hi + 1 : 0; + begin_hor_sync = fabsf(cnt_freq - 1200.0) < 50.0 ? begin_hor_sync + 1 : 0; + begin_cal_break = fabsf(cnt_freq - 1200.0) < 50.0 ? begin_cal_break + 1 : 0; + begin_cal_leader = fabsf(cal_avg - 1900.0) < 50.0 ? begin_cal_leader + 1 : 0; + begin_sep_evn = fabsf(dat_freq - 1500.0) < 50.0 ? begin_sep_evn + 1 : 0; + begin_sep_odd = fabsf(dat_freq - 2300.0) < 350.0 ? begin_sep_odd + 1 : 0; + + const float vis_tolerance = 0.9; + const float sync_tolerance = 0.7; + const float break_tolerance = 0.7; + const float leader_tolerance = 0.3; + const float seperator_tolerance = 0.7; + + int vis_ss = begin_vis_ss >= (int)(drate * vis_tolerance * vis_len) ? 1 : 0; + int vis_lo = begin_vis_lo >= (int)(drate * vis_tolerance * vis_len) ? 1 : 0; + int vis_hi = begin_vis_hi >= (int)(drate * vis_tolerance * vis_len) ? 1 : 0; + int cal_break = begin_cal_break >= (int)(drate * break_tolerance * cal_break_len) ? 1 : 0; + int cal_leader = begin_cal_leader >= (int)(drate * leader_tolerance * cal_leader_len) ? 1 : 0; + int sep_evn = begin_sep_evn >= (int)(drate * seperator_tolerance * seperator_len) ? 1 : 0; + int sep_odd = begin_sep_odd >= (int)(drate * seperator_tolerance * seperator_len) ? 1 : 0; + + // we want a pulse at the falling edge + latch_sync = begin_hor_sync > (int)(drate * sync_tolerance * hor_sync_len) ? 1 : latch_sync; + int hor_sync = begin_hor_sync > (int)(drate * sync_tolerance * hor_sync_len) ? 0 : latch_sync; + latch_sync = hor_sync ? 0 : latch_sync; + + // we only want a pulse for the bits + begin_vis_ss = vis_ss ? 0 : begin_vis_ss; + begin_vis_lo = vis_lo ? 0 : begin_vis_lo; + begin_vis_hi = vis_hi ? 0 : begin_vis_hi; + +#if DEBUG + if ((int)(ticks * dstep) < 5.0) + printf("%f %f %f %d %d %d %d %d %d %d %d\n", (float)ticks * dstep, dat_freq, cnt_freq, + 50*hor_sync+950, 50*cal_leader+850, 50*cal_break+750, + 50*vis_ss+650, 50*vis_lo+550, 50*vis_hi+450, + 50*sep_evn+350, 50*sep_odd+250); +#endif + + if (cal_leader && !cal_break && got_cal_break && + cal_ticks >= (int)(drate * (cal_leader_len + cal_break_len) * leader_tolerance) && + cal_ticks <= (int)(drate * (cal_leader_len + cal_break_len) * (2.0 - leader_tolerance))) { + vis_mode = 1; + vis_bit = -1; + dat_mode = 0; + first_hor_sync = 1; + got_cal_break = 0; + fprintf(stderr, "%f got calibration header\n", (float)ticks * dstep); + } + + if (cal_break && !cal_leader && + cal_ticks >= (int)(drate * cal_break_len * break_tolerance) && + cal_ticks <= (int)(drate * cal_break_len * (2.0 - break_tolerance))) + got_cal_break = 1; + + if (cal_leader && !cal_break) { + cal_ticks = 0; + got_cal_break = 0; + } + + if (vis_mode) { + if (vis_bit < 0) { + if (vis_ss) { + vis_ticks = 0; + vis_byte = 0; + vis_bit = 0; + dat_mode = 0; + } + } else if (vis_ticks <= (int)(drate * 10.0 * vis_len * (2.0 - vis_tolerance))) { + if (vis_ss) { + dat_mode = 1; + vis_mode = 0; + vis_bit = -1; + fprintf(stderr, "%f got VIS = 0x%x (complete)\n", (float)ticks*dstep, vis_byte); + } + if (vis_bit < 8) { + if (vis_lo) vis_bit++; + if (vis_hi) vis_byte |= 1 << vis_bit++; + } + } else { + if (vis_bit >= 8) { + dat_mode = 1; + vis_mode = 0; + vis_bit = -1; + fprintf(stderr, "%f got VIS = 0x%x (missing stop bit)\n", (float)ticks*dstep, vis_byte); + } + } + if (!vis_mode && vis_byte != 0x88) { + fprintf(stderr, "unsupported mode 0x%x, ignoring\n", vis_byte); + dat_mode = 0; + } + continue; + } + if (!dat_mode) + continue; + + // we wait until first sync + if (first_hor_sync && !hor_sync) + continue; + + // data comes after first sync + if (first_hor_sync && hor_sync) { + first_hor_sync = 0; + hor_ticks = 0; + y_pixel_x = 0; + uv_pixel_x = 0; + continue; + } + +#if DEBUG + if (hor_ticks < width) { + uint8_t *p = pixel + 3 * y * width + 3 * hor_ticks; +#if DATA + uint8_t v = limit(0.0, 255.0, 255.0 * (dat_freq - 1500.0) / 800.0); +#else + uint8_t v = limit(0.0, 255.0, 255.0 * (cnt_freq - 1100.0) / 200.0); +#endif + p[0] = v; + p[1] = v; + p[2] = v; + } +#endif + // if horizontal sync is too early, we reset to the beginning instead of ignoring + if (hor_sync && hor_ticks < (int)((hor_len - sync_porch_len) * drate)) { + for (int i = 0; i < 4; i++) { + uint8_t *p = pixel + 3 * y * width + 3 * (width - i - 10); + p[0] = 255; + p[1] = 0; + p[2] = 255; + } + hor_ticks = 0; + y_pixel_x = 0; + uv_pixel_x = 0; + } + + // we always sync if sync pulse is where it should be. + if (hor_sync && (hor_ticks >= (int)((hor_len - sync_porch_len) * drate) && + hor_ticks < (int)((hor_len + sync_porch_len) * drate))) { +#if DEBUG + uint8_t *p = pixel + 3 * y * width + 3 * hor_ticks + 6 * (int)(sync_porch_len * drate); + p[0] = 0; + p[1] = 255; + p[2] = 255; + y++; +#else + process_line(pixel, y_pixel, uv_pixel, y_width, uv_width, width, height, y++); +#endif + if (y == height) + break; + odd ^= 1; + hor_ticks = 0; + y_pixel_x = 0; + uv_pixel_x = 0; + } + + // if horizontal sync is missing, we extrapolate from last sync + if (hor_ticks >= (int)((hor_len + sync_porch_len) * drate)) { +#if DEBUG + for (int i = 0; i < 4; i++) { + uint8_t *p = pixel + 3 * y * width + 3 * (width - i - 5); + p[0] = 255; + p[1] = 255; + p[2] = 0; + } + y++; +#else + process_line(pixel, y_pixel, uv_pixel, y_width, uv_width, width, height, y++); +#endif + if (y == height) + break; + odd ^= 1; + missing_sync++; + hor_ticks -= (int)(hor_len * drate); + // we are not at the pixels yet, so no correction here + y_pixel_x = 0; + uv_pixel_x = 0; + } + + if (hor_ticks > (int)((sync_porch_len + y_len) * drate) && hor_ticks < (int)((sync_porch_len + y_len + seperator_len) * drate)) { + odd_count += sep_odd; + evn_count += sep_evn; + } + // we try to correct from odd / even seperator + if (evn_count != odd_count && hor_ticks > (int)((sync_porch_len + y_len + seperator_len) * drate)) { + // even seperator + if (evn_count > odd_count && odd) { + odd = 0; + seperator_correction++; +#if DEBUG + for (int i = 0; i < 4; i++) { + uint8_t *p = pixel + 3 * y * width + 3 * (width - i - 15); + p[0] = 255; + p[1] = 0; + p[2] = 0; + } +#endif + } + // odd seperator + if (odd_count > evn_count && !odd) { + odd = 1; + seperator_correction++; +#if DEBUG + for (int i = 0; i < 4; i++) { + uint8_t *p = pixel + 3 * y * width + 3 * (width - i - 15); + p[0] = 0; + p[1] = 255; + p[2] = 0; + } +#endif + } + evn_count = 0; + odd_count = 0; + } +#if DEBUG + float fixme = 0.0007; + if (hor_ticks == (int)((fixme + sync_porch_len) * drate) || + hor_ticks == (int)((fixme + sync_porch_len + y_len) * drate) || + hor_ticks == (int)((fixme + sync_porch_len + y_len + seperator_len) * drate) || + hor_ticks == (int)((fixme + sync_porch_len + y_len + seperator_len + porch_len) * drate) || + hor_ticks == (int)((fixme + sync_porch_len + y_len + seperator_len + porch_len + uv_len) * drate)) { + uint8_t *p = pixel + 3 * y * width + 3 * hor_ticks; + p[0] = 255; + p[1] = 0; + p[2] = 0; + } +#else + // TODO: need better way to compensate for pulse decay time + float fixme = 0.0007; + if (y_pixel_x < y_width && hor_ticks >= (int)((fixme + sync_porch_len) * drate)) + y_pixel[y_pixel_x++ + (y % 2) * y_width] = limit(0.0, 255.0, 255.0 * (dat_freq - 1500.0) / 800.0); + + if (uv_pixel_x < uv_width && hor_ticks >= (int)((fixme + sync_porch_len + y_len + seperator_len + porch_len) * drate)) + uv_pixel[uv_pixel_x++ + odd * uv_width] = limit(0.0, 255.0, 255.0 * (dat_freq - 1500.0) / 800.0); +#endif + } + + munmap_file(wav_p, wav_size); + + free_ddc(cnt_ddc); + free_ddc(dat_ddc); + free(cnt_amp); + free(dat_amp); + + munmap_file(ppm_p, ppm_size); + fprintf(stderr, "%d missing sync's and %d corrections from seperator\n", missing_sync, seperator_correction); + + return 0; +} + diff --git a/encode.c b/encode.c new file mode 100644 index 0000000..3c349b8 --- /dev/null +++ b/encode.c @@ -0,0 +1,355 @@ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +float limit(float min, float max, float x) +{ + float tmp = x < min ? min : x; + return tmp > max ? max : tmp; +} +float lerp(float a, float b, float x) +{ + return a - a * x + b * x; +} +uint8_t R_YUV(uint8_t Y, uint8_t U, uint8_t V) +{ + (void)U; + return limit(0.0, 255.0, 0.003906 * ((298.082 * (Y - 16.0)) + (408.583 * (V - 128.0)))); +} +uint8_t G_YUV(uint8_t Y, uint8_t U, uint8_t V) +{ + return limit(0.0, 255.0, 0.003906 * ((298.082 * (Y - 16.0)) + (-100.291 * (U - 128.0)) + (-208.12 * (V - 128.0)))); +} +uint8_t B_YUV(uint8_t Y, uint8_t U, uint8_t V) +{ + (void)V; + return limit(0.0, 255.0, 0.003906 * ((298.082 * (Y - 16.0)) + (516.411 * (U - 128.0)))); +} + +uint8_t Y_RGB(uint8_t R, uint8_t G, uint8_t B) +{ + return limit(0.0, 255.0, 16.0 + (0.003906 * ((65.738 * R) + (129.057 * G) + (25.064 * B)))); +} +uint8_t V_RGB(uint8_t R, uint8_t G, uint8_t B) +{ + return limit(0.0, 255.0, 128.0 + (0.003906 * ((112.439 * R) + (-94.154 * G) + (-18.285 * B)))); +} +uint8_t U_RGB(uint8_t R, uint8_t G, uint8_t B) +{ + return limit(0.0, 255.0, 128.0 + (0.003906 * ((-37.945 * R) + (-74.494 * G) + (112.439 * B)))); +} + +void *mmap_file_ro(char *name, size_t *size) +{ + *size = 0; + int fd = open(name, O_RDONLY); + if (fd == -1) { + perror("open"); + return 0; + } + + struct stat sb; + if (fstat(fd, &sb) == -1) { + perror("fstat"); + return 0; + } + + if (!S_ISREG(sb.st_mode)) { + fprintf(stderr, "%s not a file\n", name); + return 0; + } + + void *p = mmap(0, sb.st_size, PROT_READ, MAP_SHARED, fd, 0); + if (p == MAP_FAILED) { + perror("mmap"); + return 0; + } + + if (close(fd) == -1) { + perror ("close"); + return 0; + } + *size = sb.st_size; + return p; +} + +void *mmap_file_rw(char *name, size_t size) +{ + int fd = open(name, O_RDWR|O_CREAT|O_TRUNC, S_IRUSR|S_IWUSR|S_IRGRP|S_IROTH); + if (fd == -1) { + perror("open"); + return 0; + } + + struct stat sb; + if (fstat(fd, &sb) == -1) { + perror("fstat"); + return 0; + } + + if (!S_ISREG(sb.st_mode)) { + fprintf(stderr, "%s not a file\n", name); + return 0; + } + + if (lseek(fd, size - 1, SEEK_SET) == -1) { + perror("lseek"); + return 0; + } + + if (write(fd, "", 1) != 1) { + perror("write"); + return 0; + } + + void *p = mmap(0, size, PROT_READ|PROT_WRITE, MAP_SHARED, fd, 0); + if (p == MAP_FAILED) { + perror("mmap"); + return 0; + } + + if (close(fd) == -1) { + perror ("close"); + return 0; + } + return p; +} + +int munmap_file(void *p, size_t size) +{ + if (munmap(p, size) == -1) { + perror("munmap"); + return 0; + } + return 1; +} + +typedef struct { + uint32_t ChunkID; + uint32_t ChunkSize; + uint32_t Format; + uint32_t Subchunk1ID; + uint32_t Subchunk1Size; + uint16_t AudioFormat; + uint16_t NumChannels; + uint32_t SampleRate; + uint32_t ByteRate; + uint16_t BlockAlign; + uint16_t BitsPerSample; + uint32_t Subchunk2ID; + uint32_t Subchunk2Size; +} wav_t; + +short *buffer; +complex float nco; +int sample; +float hz2rad; + +void add_sample(float val) { +// static float avg = 0.0; +// const float a = 0.9; +// avg = a * val + (1.0 - a) * avg; +// buffer[sample++] = (float)SHRT_MAX * avg; + buffer[sample++] = (float)SHRT_MAX * val; +// buffer[sample++] = (float)SHRT_MAX * avg + random() / (RAND_MAX / 10000); +} +void add_freq(float freq) { + add_sample(creal(nco)); + nco *= cexpf(freq * hz2rad * I); +} + +int main(int argc, char **argv) +{ + if (argc != 4) { + fprintf(stderr, "usage: %s \n", argv[0]); + return 1; + } + + size_t ppm_size; + char *ppm_p = mmap_file_ro(argv[1], &ppm_size); + const int width = 320; + const int height = 240; + const char *ppm_head = "P6 320 240 255\n"; + + if (strncmp(ppm_head, ppm_p, strlen(ppm_head))) { + fprintf(stderr, "unsupported image file\n"); + return 1; + } + + uint8_t *pixel = (uint8_t *)ppm_p + strlen(ppm_head); + + float rate = atoi(argv[3]); + + if (fabsf(0.0015 * rate - (int)(0.0015 * rate)) > 0.0001) + fprintf(stderr, "this rate will not give accurate (smooth) results.\ntry 40000Hz and resample to %0.fHz\n", rate); + + hz2rad = (2.0 * M_PI) / rate; + nco = -I * 0.7; + enum { N = 13 }; + float seq_freq[N] = { 1900.0, 1200.0, 1900.0, 1200.0, 1300.0, 1300.0, 1300.0, 1100.0, 1300.0, 1300.0, 1300.0, 1100.0, 1200.0 }; + float seq_time[N] = { 0.3, 0.01, 0.3, 0.03, 0.03, 0.03, 0.03, 0.03, 0.03, 0.03, 0.03, 0.03, 0.03 }; + + size_t wav_size = 4096 * ((size_t)(37.5 * rate * 2 + 44 + 4095) / 4096); + int samples = (wav_size - 44) / 2; + char *wav_p = mmap_file_rw(argv[2], wav_size); + if (!wav_p) + return 1; + + buffer = (short *)(wav_p + sizeof(wav_t)); + + sample = 0; + + for (int ticks = 0; ticks < (int)(0.3 * rate); ticks++) + add_sample(0.0); + + for (int i = 0; i < N; i++) + for (int ticks = 0; ticks < (int)(seq_time[i] * rate); ticks++) + add_freq(seq_freq[i]); + + for (int y = 0; y < height; y++) { + // EVEN LINES + // SYNC + for (int ticks = 0; ticks < (int)(0.009 * rate); ticks++) { + add_freq(1200.0); + } + // PORCH + for (int ticks = 0; ticks < (int)(0.003 * rate); ticks++) { + add_freq(1500.0); + } + // Y + for (int ticks = 0; ticks < (int)(0.088 * rate); ticks++) { + float xf = limit(0.0, 319.0, (320.0 * (float)ticks) / (0.088 * rate)); + int x0 = xf; + int x1 = limit(0.0, 319.0, x0 + 1); + int off0 = 3 * y * width + 3 * x0; + int off1 = 3 * y * width + 3 * x1; + uint8_t R0 = pixel[off0 + 0]; + uint8_t G0 = pixel[off0 + 1]; + uint8_t B0 = pixel[off0 + 2]; + uint8_t R1 = pixel[off1 + 0]; + uint8_t G1 = pixel[off1 + 1]; + uint8_t B1 = pixel[off1 + 2]; + uint8_t R = lerp(R0, R1, xf - (float)x0); + uint8_t G = lerp(G0, G1, xf - (float)x0); + uint8_t B = lerp(B0, B1, xf - (float)x0); + add_freq(1500.0 + 800.0 * Y_RGB(R, G, B) / 255.0); + } + // EVEN + for (int ticks = 0; ticks < (int)(0.0045 * rate); ticks++) { + add_freq(1500.0); + } + // PORCH + for (int ticks = 0; ticks < (int)(0.0015 * rate); ticks++) { + add_freq(1900.0); + } + // V + for (int ticks = 0; ticks < (int)(0.044 * rate); ticks++) { + float xf = limit(0.0, 159.0, (160.0 * (float)ticks) / (0.044 * rate)); + int x0 = xf; + int x1 = limit(0.0, 159.0, x0 + 1); + int evn0 = 3 * y * width + 6 * x0; + int evn1 = 3 * y * width + 6 * x1; + int odd0 = 3 * (y + 1) * width + 6 * x0; + int odd1 = 3 * (y + 1) * width + 6 * x1; + uint8_t R0 = (pixel[evn0 + 0] + pixel[odd0 + 0] + pixel[evn0 + 3] + pixel[odd0 + 3]) / 4; + uint8_t G0 = (pixel[evn0 + 1] + pixel[odd0 + 1] + pixel[evn0 + 4] + pixel[odd0 + 4]) / 4; + uint8_t B0 = (pixel[evn0 + 2] + pixel[odd0 + 2] + pixel[evn0 + 5] + pixel[odd0 + 5]) / 4; + uint8_t R1 = (pixel[evn1 + 0] + pixel[odd1 + 0] + pixel[evn1 + 3] + pixel[odd1 + 3]) / 4; + uint8_t G1 = (pixel[evn1 + 1] + pixel[odd1 + 1] + pixel[evn1 + 4] + pixel[odd1 + 4]) / 4; + uint8_t B1 = (pixel[evn1 + 2] + pixel[odd1 + 2] + pixel[evn1 + 5] + pixel[odd1 + 5]) / 4; + uint8_t R = lerp(R0, R1, xf - (float)x0); + uint8_t G = lerp(G0, G1, xf - (float)x0); + uint8_t B = lerp(B0, B1, xf - (float)x0); + add_freq(1500.0 + 800.0 * V_RGB(R, G, B) / 255.0); + } + // ODD LINES + y++; + // SYNC + for (int ticks = 0; ticks < (int)(0.009 * rate); ticks++) { + add_freq(1200.0); + } + // PORCH + for (int ticks = 0; ticks < (int)(0.003 * rate); ticks++) { + add_freq(1500.0); + } + // Y + for (int ticks = 0; ticks < (int)(0.088 * rate); ticks++) { + float xf = limit(0.0, 319.0, (320.0 * (float)ticks) / (0.088 * rate)); + int x0 = xf; + int x1 = limit(0.0, 319.0, x0 + 1); + int off0 = 3 * y * width + 3 * x0; + int off1 = 3 * y * width + 3 * x1; + uint8_t R0 = pixel[off0 + 0]; + uint8_t G0 = pixel[off0 + 1]; + uint8_t B0 = pixel[off0 + 2]; + uint8_t R1 = pixel[off1 + 0]; + uint8_t G1 = pixel[off1 + 1]; + uint8_t B1 = pixel[off1 + 2]; + uint8_t R = lerp(R0, R1, xf - (float)x0); + uint8_t G = lerp(G0, G1, xf - (float)x0); + uint8_t B = lerp(B0, B1, xf - (float)x0); + add_freq(1500.0 + 800.0 * Y_RGB(R, G, B) / 255.0); + } + // ODD + for (int ticks = 0; ticks < (int)(0.0045 * rate); ticks++) { + add_freq(2300.0); + } + // PORCH + for (int ticks = 0; ticks < (int)(0.0015 * rate); ticks++) { + add_freq(1900.0); + } + // U + for (int ticks = 0; ticks < (int)(0.044 * rate); ticks++) { + float xf = limit(0.0, 159.0, (160.0 * (float)ticks) / (0.044 * rate)); + int x0 = xf; + int x1 = limit(0.0, 159.0, x0 + 1); + int evn0 = 3 * (y - 1) * width + 6 * x0; + int evn1 = 3 * (y - 1) * width + 6 * x1; + int odd0 = 3 * y * width + 6 * x0; + int odd1 = 3 * y * width + 6 * x1; + uint8_t R0 = (pixel[evn0 + 0] + pixel[odd0 + 0] + pixel[evn0 + 3] + pixel[odd0 + 3]) / 4; + uint8_t G0 = (pixel[evn0 + 1] + pixel[odd0 + 1] + pixel[evn0 + 4] + pixel[odd0 + 4]) / 4; + uint8_t B0 = (pixel[evn0 + 2] + pixel[odd0 + 2] + pixel[evn0 + 5] + pixel[odd0 + 5]) / 4; + uint8_t R1 = (pixel[evn1 + 0] + pixel[odd1 + 0] + pixel[evn1 + 3] + pixel[odd1 + 3]) / 4; + uint8_t G1 = (pixel[evn1 + 1] + pixel[odd1 + 1] + pixel[evn1 + 4] + pixel[odd1 + 4]) / 4; + uint8_t B1 = (pixel[evn1 + 2] + pixel[odd1 + 2] + pixel[evn1 + 5] + pixel[odd1 + 5]) / 4; + uint8_t R = lerp(R0, R1, xf - (float)x0); + uint8_t G = lerp(G0, G1, xf - (float)x0); + uint8_t B = lerp(B0, B1, xf - (float)x0); + add_freq(1500.0 + 800.0 * U_RGB(R, G, B) / 255.0); + } + } + + while (sample < samples) + add_sample(0.0); + + wav_t *wav = (wav_t *)wav_p; + wav->ChunkID = 0x46464952; + wav->ChunkSize = 36 + 2 * samples; + wav->Format = 0x45564157; + wav->Subchunk1ID = 0x20746d66; + wav->Subchunk1Size = 16; + wav->AudioFormat = 1; + wav->NumChannels = 1; + wav->SampleRate = rate; + wav->ByteRate = 2 * rate; + wav->BlockAlign = 2; + wav->BitsPerSample = 16; + wav->Subchunk2ID = 0x61746164; + wav->Subchunk2Size = 2 * samples; + + munmap_file(wav_p, wav_size); + munmap_file(ppm_p, ppm_size); + return 0; +} + diff --git a/smpte.ppm b/smpte.ppm new file mode 100644 index 0000000000000000000000000000000000000000..5cd63d97fd2615ccc3b633f6e5d2dda20caae020 GIT binary patch literal 230416 zcmeI%zlvT}7=`gcjO0%dArP>T#LgmxL?c*O*;;rLTG$G68OF-aEAS@18&NMsPfWA3 zR*)p@?b!pr!w6?BbkBI^d(N}pJ>Px#`17ZqJ%0Md@z1kocQ1Qfzy0>Z<#M&d*QdY! zw8hVHxxDxL2UnM?)4v}5@@S9W%;ncx-)!(aUi?BHy#Mq0znFhN8GGP+!Ee7b@Eh+9 z{Ps%&zn51(CHx!4=iin-|HgZRfBE;-_h`IUUqt@p-`9LBUj5wgFaNf1qv2ovZQ({^ z{(be$%D?dzc@a0c?4JYPBjOjo2ciG{h2Pjw+ds!g#`#qvWd9uR4gTfdux&@@SMB29 z=&Xoe{5G)B;1|EKTkdH8Z2#OF=4kAnL-hR1zhT>s{L8;B+-Ufhf5WyN`ImoNxY79g zQ7D*y`8RCak$?HOg&Ph3@^9F-BmeSm3pbkJ=&b1c3co>XNBif1_lTQZ_{DDn8x4N( z+XHVj=3nz~OuHS;zu|rU<=?PvNB-sCux&^F<=+-=H2llIE!=3FKZk<(mw#Kh(eN+- zhHX3YFaL&ZJMu68ws4~fj?RkCukahRcC>#Ec#pWrgOr0~-x~@!JD$H0EFPZ%n%#&A;J&{^j4WZAbp) z->_{*{^j2mZZ!PMzb)KooIi(x`ImoNxY6)0|AuWl@-P2}Z9DQW|F&?W369Q+&adzr zw05+A4tS5a$%SA1Hn7p)7r!y>b~OK*e|zAK#{3)J=U@H}+jit%{%zq#!@v9+w(ZEj z{M*8f#`$w7n1A^N-sfNb4cm6)U;b_3M#I1S8@BDpzx><6jmG(ND42iwH*DLHfBCnC8x8;R zZ`igY|MG7OH=5w+tmym-zd>t9`{#i7h?`va#cu-}4Sw<418+3uU-NHFyB*EH;eGz) z->_{*{^j4WZAbp)-xh8({L8;B+-RIXhl2T+e_OcG@Gt*{Z9DQW|AuWl@-P3kaH9#1 z&Wg^j@Ef#tw0{nGkGRQ&U;H+((cl-qG3|CV|C)b$;El%o8{X$%{ter9XNBif1 z_lTQZ_{DDn8x4N(+XHVj=3nz~OuHS;zu|rU<=?PvNB-sCux&^F<=+-=H2llIE!=3F zKZk<(mw#Kh(eN+-hHX3YFaL&ZJMu68ws4~fj?RkCukahRcC>#Ec#pWrgOr0~-x~@!JD$H0EFPZ%n%# z&A;J&{^j4WZAbp)->_{*{^j2mZZ!PMzb)KooIi(x`ImoNxY6)0|AuWl@-P2}Z9DQW z|F&?W369Q+&adzrw05+A4tS5a$%SA1Hn7p)7r!y>b~OK*e|zAK#{3)J=U@H}+jit% z{%zq#!@v9+w(ZEj{M*8f#`$w7n1A^N-sfNb4cm6)U;b_3M#I1S8@BDpzx><6jmG(ND42iw zH*DLHfBCnC8x8;RZ`igY|MG7OH=5w+tmym-zd>t9`{#i7h?`va#cu-}4Sw<418+3u zU-NHFyB*EH;eGz)->_{*{^j4WZAbp)-xh8({L8;B+-RIXhl2T+e_OcG@Gt*{Z9DQW z|AuWl@-P3kaH9#1&Wg^j@Ef#tw0{nGkGRQ&U;H+((cl-qG3|CV|C)b$;El%o8{X$% z{ter9XNBif1_lTQZ_{DDn8x4N(+XHVj=3nz~OuHS;zu|rU<=?PvNB-sCux&^F z<=+-=H2llIE!=3FKZk<(mw#Kh(eN+-hHX3YFaL&ZJMu68ws4~fj?RkCukahRcC>#E zc#pWrg37R1&F$;!9UjQ>GdzF({IbWj3C{TKw*-FUErH*DXW%zp>HHfHgn#3K@Nd5~_?Lgf zwjKGGfBRh;4gbai;oo>5M>PM2g87$!Te#8iFaP$tG#dVm2g1McK#s`2{M*8f#`!b< zws4~fj?Rkq&jIfd@r&P}wc}a9#`nc9euLJI=3o5wEEA3SHy()jHy+3l&A-?dIU@h^Z`igY|MG7OH=5w+tmym-zd>t9`{#i7h-U#C-xt65 zZD6CpFMfNLiN^dJ55)W%59EmEU;b_3Mq~cv-xh8(!O>aK{yE@1B7X53w01lT*!aHq z#c$Bs(fs=wen0v4>V|y%@#kNE-;nF;kG_9;_2Z5B`sChMJ>Tuw+`oVS&KqC%?%jK< z|BH{`=6d+>;ae@f?q9ws@E$)WepkE@@cR<3x54*g{Jz6`1HbrvWBuOY-*-T4>3@P> z{NgvZ;!BB}yEuOFi(mZy4`*P`ui|4iB|5+A$2eD!OIPDx{^j35cqwsn7soGt@rz&k z=kYz@??*3*@b{y+w(ZsUcjRu?J$~`qtgUG9JB!(I9sJ^V1+Bd%en;+R-QyR(&Dx3v zzq6Pf*TFA-SJ2vP;&myaUJ~PcLlAzCVofmX5Hf#zs=f;2EVhI z9oNAxepk@iYvOn0Zq_}1@!PDeXz)9W*>N5G;&%nDy(WG~?q=QN7r)KgiUz;4m>t)_ zFMe0h+H2x>myaUJ~PcLlAzCVofmX5Hf#zs=f;2EVhI9oNAxepk@iYvOn0 zZq_}1@!PDeXz)9W*>N5G;&%nDy(WG~?q=QN7r)KgiUz;4m>t)_FMe0h+H2x>myaUJ~PcLlAzCVofmX5Hf#zs=f;2EVhI9oNAxepk@iYvOn0Zq_}1@!PDeXz)9W z*>N5G;&%nDy(WG~?q=QN7r)KgiUz;4m>t)_FMe0h+H2x>