mirror
/
slowrx
kopia lustrzana https://github.com/windytan/slowrx
1
0
Forkuj 0
Kod Zgłoszenia Projekty Wydania Wiki Aktywność

better frequency estimator using complex-Hermitian DFT

pull/5/head
Oona 2013-01-26 18:38:40 +02:00
rodzic 6369f9cd5e
commit 9132c32227
8 zmienionych plików z 57 dodań i 53 usunięć
Pokaż statystyki Ściągnij plik aktualizacji Ściągnij plik porównania Expand all files Collapse all files

8
common.c
Wyświetl plik

@ -13,16 +13,15 @@
gboolean Abort = FALSE;
gboolean Adaptive = TRUE;
double *in = NULL;
gboolean *HasSync = NULL;
gshort HedrShift = 0;
gboolean ManualActivated = FALSE;
gboolean ManualResync = FALSE;
double *out = NULL;
guchar *StoredLum = NULL;
pthread_t thread1;
FFTStuff fft;
GuiObjs gui;
PicMeta CurrentPic;
PcmData pcm;
@ -44,6 +43,11 @@ guint GetBin (double Freq, guint FFTLen) {
return (Freq / 44100 * FFTLen);
}
// Sinusoid power from complex DFT coefficients
double power (fftw_complex coeff) {
return pow(coeff[0],2) + pow(coeff[1],2);
}
// Clip to [0..255]
guchar clip (double a) {
if (a < 0) return 0;

11
common.h
Wyświetl plik

@ -9,18 +9,22 @@
extern gboolean Abort;
extern gboolean Adaptive;
extern gboolean *HasSync;
extern double *in;
extern gboolean ManualActivated;
extern gboolean ManualResync;
extern double *out;
extern guchar *StoredLum;
extern pthread_t thread1;
extern guchar VISmap[];
typedef struct _FFTStuff FFTStuff;
struct _FFTStuff {
double *in;
fftw_complex *out;
};
extern FFTStuff fft;
typedef struct _PcmData PcmData;
struct _PcmData {
snd_pcm_t *handle;
int PeakVal;
gint16 *Buffer;
int WindowPtr;
gboolean BufferDrop;
@ -119,6 +123,7 @@ typedef struct ModeSpec {
extern _ModeSpec ModeSpec[];
double power (fftw_complex coeff);
guchar clip (double a);
void createGUI ();
double deg2rad (double Deg);

8
fsk.c
Wyświetl plik

@ -35,7 +35,7 @@ void GetFSK (char *dest) {
0x03, 0x23, 0x13, 0x33, 0x0b, 0x2b, 0x1b, 0x3b,
0x07, 0x27, 0x17, 0x37, 0x0f, 0x2f, 0x1f, 0x3f };
for (i = 0; i < FFTLen; i++) in[i] = 0;
for (i = 0; i < FFTLen; i++) fft.in[i] = 0;
// Create 22ms Hann window
for (i = 0; i < 970; i++) Hann[i] = 0.5 * (1 - cos( 2 * M_PI * i / 969.0 ) );
@ -46,7 +46,7 @@ void GetFSK (char *dest) {
readPcm(InSync ? 970: 485);
// Apply Hann window
for (i = 0; i < 970; i++) in[i] = pcm.Buffer[pcm.WindowPtr+i- 485] * Hann[i];
for (i = 0; i < 970; i++) fft.in[i] = pcm.Buffer[pcm.WindowPtr+i- 485] * Hann[i];
pcm.WindowPtr += (InSync ? 970 : 485);
@ -61,8 +61,8 @@ void GetFSK (char *dest) {
HiPow = 0;
for (i = LoBin; i <= HiBin; i++) {
if (i < MidBin) LoPow += pow(out[i], 2) + pow(out[FFTLen - i], 2);
else HiPow += pow(out[i], 2) + pow(out[FFTLen - i], 2);
if (i < MidBin) LoPow += power(fft.out[i]);
else HiPow += power(fft.out[i]);
}
Bit = (LoPow>HiPow);

1
gui.c
Wyświetl plik

@ -86,7 +86,6 @@ void createGUI() {
// Draw signal level meters according to given values
void setVU (double *Power, int FFTLen, int WinIdx) {
int x,y;
int PWRdB;
guchar *pixelsPWR, *pixelsSNR, *pPWR, *pSNR;
unsigned int rowstridePWR,rowstrideSNR;
double logpow;

12
pcm.c
Wyświetl plik

@ -53,12 +53,10 @@ void readPcm(gint numsamples) {
pcm.Buffer[i] = tmp[i] & 0xffff;
pcm.WindowPtr = BUFLEN/2;
} else {
for (i=0; i<BUFLEN-numsamples; i++) pcm.Buffer[i] = pcm.Buffer[i+numsamples];
for (i=0; i<numsamples; i++) {
pcm.Buffer[BUFLEN-numsamples+i] = tmp[i] & 0xffff;
// Keep track of max power for VU meter
if (abs(pcm.Buffer[i]) > pcm.PeakVal) pcm.PeakVal = abs(pcm.Buffer[i]);
}
// Move buffer and push samples
for (i=0; i<BUFLEN-numsamples; i++) pcm.Buffer[i] = pcm.Buffer[i+numsamples];
for (i=BUFLEN-numsamples; i<BUFLEN; i++) pcm.Buffer[i] = tmp[i-(BUFLEN-numsamples)] & 0xffff;
pcm.WindowPtr -= numsamples;
}
@ -110,7 +108,7 @@ int initPcmDevice(char *wanteddevname) {
char pcm_name[30];
unsigned int exact_rate = 44100;
int card;
gboolean found;
gboolean found;
char *cardname;
pcm.BufferDrop = FALSE;

23
slowrx.c
Wyświetl plik

@ -28,7 +28,7 @@ void *Listen() {
guchar Mode=0;
struct tm *timeptr = NULL;
time_t timet;
gboolean Finished;
gboolean Finished;
char id[20];
GtkTreeIter iter;
@ -215,22 +215,21 @@ int main(int argc, char *argv[]) {
}
// Prepare FFT
in = fftw_malloc(sizeof(double) * 2048);
if (in == NULL) {
fft.in = fftw_alloc_real(2048);
if (fft.in == NULL) {
perror("GetVideo: Unable to allocate memory for FFT");
exit(EXIT_FAILURE);
}
out = fftw_malloc(sizeof(double) * 2048);
if (out == NULL) {
fft.out = fftw_alloc_complex(2048);
if (fft.out == NULL) {
perror("GetVideo: Unable to allocate memory for FFT");
fftw_free(in);
fftw_free(fft.in);
exit(EXIT_FAILURE);
}
memset(in, 0, sizeof(double) * 2048);
memset(out, 0, sizeof(double) * 2048);
memset(fft.in, 0, sizeof(double) * 2048);
Plan1024 = fftw_plan_r2r_1d(1024, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
Plan2048 = fftw_plan_r2r_1d(2048, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
Plan1024 = fftw_plan_dft_r2c_1d(1024, fft.in, fft.out, FFTW_ESTIMATE);
Plan2048 = fftw_plan_dft_r2c_1d(2048, fft.in, fft.out, FFTW_ESTIMATE);
createGUI();
populateDeviceList();
@ -250,8 +249,8 @@ int main(int argc, char *argv[]) {
g_object_unref(pixbuf_rx);
free(StoredLum);
fftw_free(in);
fftw_free(out);
fftw_free(fft.in);
fftw_free(fft.out);
return (EXIT_SUCCESS);
}

40
video.c
Wyświetl plik

@ -54,8 +54,10 @@ gboolean GetVideo(guchar Mode, double Rate, int Skip, gboolean Redraw) {
gushort HannLens[7] = { 48, 64, 96, 128, 256, 512, 1024 };
for (j = 0; j < 7; j++)
for (i = 0; i < HannLens[j]; i++)
//Hann[j][i] = exp(-0.5*pow((i-(HannLens[j]-1)*0.5)/(0.35*(HannLens[j]-1)*0.5),2));
Hann[j][i] = 0.5 * (1 - cos( (2 * M_PI * i) / (HannLens[j] - 1)) );
// Starting times of video channels on every line, counted from beginning of line
switch (Mode) {
@ -194,20 +196,19 @@ gboolean GetVideo(guchar Mode, double Rate, int Skip, gboolean Redraw) {
Praw = Psync = 0;
memset(in, 0, sizeof(in[0]) *FFTLen);
memset(out, 0, sizeof(out[0])*FFTLen);
memset(fft.in, 0, sizeof(double)*FFTLen);
// Hann window
for (i = 0; i < 64; i++) in[i] = pcm.Buffer[pcm.WindowPtr+i-32] / 32768.0 * Hann[1][i];
for (i = 0; i < 64; i++) fft.in[i] = pcm.Buffer[pcm.WindowPtr+i-32] / 32768.0 * Hann[1][i];
fftw_execute(Plan1024);
for (i=0;i<LopassBin;i++) {
if (i >= GetBin(1500+CurrentPic.HedrShift, FFTLen) && i <= GetBin(2300+CurrentPic.HedrShift, FFTLen))
Praw += pow(out[i], 2) + pow(out[FFTLen-i], 2);
Praw += power(fft.out[i]);
if (i >= SyncTargetBin-1 && i <= SyncTargetBin+1)
Psync += (pow(out[i], 2) + pow(out[FFTLen-i], 2)) * (1- .5*abs(SyncTargetBin-i));
Psync += power(fft.out[i]) * (1- .5*abs(SyncTargetBin-i));
}
Praw /= (GetBin(2300+CurrentPic.HedrShift, FFTLen) - GetBin(1500+CurrentPic.HedrShift, FFTLen));
@ -229,29 +230,28 @@ gboolean GetVideo(guchar Mode, double Rate, int Skip, gboolean Redraw) {
if (SampleNum == NextSNRtime) {
memset(in, 0, sizeof(double)*FFTLen);
memset(out, 0, sizeof(double)*FFTLen);
memset(fft.in, 0, sizeof(double)*FFTLen);
memset(fft.out, 0, sizeof(double)*FFTLen);
// Apply Hann window
for (i = 0; i < FFTLen; i++) in[i] = pcm.Buffer[pcm.WindowPtr + i - FFTLen/2] / 32768.0 * Hann[6][i];
for (i = 0; i < FFTLen; i++) fft.in[i] = pcm.Buffer[pcm.WindowPtr + i - FFTLen/2] / 32768.0 * Hann[6][i];
// FFT
fftw_execute(Plan1024);
// Calculate video-plus-noise power (1500-2300 Hz)
Pvideo_plus_noise = 0;
for (n = GetBin(1500+CurrentPic.HedrShift, FFTLen); n <= GetBin(2300+CurrentPic.HedrShift, FFTLen); n++)
Pvideo_plus_noise += pow(out[n], 2) + pow(out[FFTLen - n], 2);
Pvideo_plus_noise += power(fft.out[n]);
// Calculate noise-only power (400-800 Hz + 2700-3400 Hz)
Pnoise_only = 0;
for (n = GetBin(400+CurrentPic.HedrShift, FFTLen); n <= GetBin(800+CurrentPic.HedrShift, FFTLen); n++)
Pnoise_only += pow(out[n], 2) + pow(out[FFTLen - n], 2);
Pnoise_only += power(fft.out[n]);
for (n = GetBin(2700+CurrentPic.HedrShift, FFTLen); n <= GetBin(3400+CurrentPic.HedrShift, FFTLen); n++)
Pnoise_only += pow(out[n], 2) + pow(out[FFTLen - n], 2);
Pnoise_only += power(fft.out[n]);
// Bandwidths
VideoPlusNoiseBins = GetBin(2300, FFTLen) - GetBin(1500, FFTLen) + 1;
@ -294,14 +294,13 @@ gboolean GetVideo(guchar Mode, double Rate, int Skip, gboolean Redraw) {
// Minimum winlength can be doubled for Scottie DX
if (Mode == SDX && WinIdx < 6) WinIdx++;
WinLength = HannLens[WinIdx];
memset(in, 0, sizeof(double)*FFTLen);
memset(out, 0, sizeof(double)*FFTLen);
memset(fft.in, 0, sizeof(double)*FFTLen);
memset(Power, 0, sizeof(double)*1024);
// Apply window function
for (i = 0; i < WinLength; i++) in[i] = pcm.Buffer[pcm.WindowPtr + i - WinLength/2] / 32768.0 * Hann[WinIdx][i];
WinLength = HannLens[WinIdx];
for (i = 0; i < WinLength; i++) fft.in[i] = pcm.Buffer[pcm.WindowPtr + i - WinLength/2] / 32768.0 * Hann[WinIdx][i];
fftw_execute(Plan1024);
@ -310,7 +309,7 @@ gboolean GetVideo(guchar Mode, double Rate, int Skip, gboolean Redraw) {
// Find the bin with most power
for (n = GetBin(1500 + CurrentPic.HedrShift, FFTLen) - 1; n <= GetBin(2300 + CurrentPic.HedrShift, FFTLen) + 1; n++) {
Power[n] = pow(out[n],2) + pow(out[FFTLen - n], 2);
Power[n] = power(fft.out[n]);
if (MaxBin == 0 || Power[n] > Power[MaxBin]) MaxBin = n;
}
@ -326,6 +325,8 @@ gboolean GetVideo(guchar Mode, double Rate, int Skip, gboolean Redraw) {
Freq = ( (MaxBin > GetBin(1900 + CurrentPic.HedrShift, FFTLen)) ? 2300 : 1500 ) + CurrentPic.HedrShift;
}
printf("%.0f\n",Freq);
} /* endif (SampleNum == PixelGrid[PixelIdx].Time) */
@ -404,7 +405,6 @@ gboolean GetVideo(guchar Mode, double Rate, int Skip, gboolean Redraw) {
if (!Redraw && SampleNum % 8820 == 0) {
setVU(Power, FFTLen, WinIdx);
pcm.PeakVal = 0;
}
if (Abort) return FALSE;

7
vis.c
Wyświetl plik

@ -24,7 +24,7 @@ guchar GetVIS () {
gboolean gotvis = FALSE;
guchar Bit[8] = {0}, ParityBit = 0;
for (i = 0; i < FFTLen; i++) in[i] = 0;
for (i = 0; i < FFTLen; i++) fft.in[i] = 0;
// Create 20ms Hann window
for (i = 0; i < 882; i++) Hann[i] = 0.5 * (1 - cos( (2 * M_PI * (double)i) / 881 ) );
@ -45,7 +45,7 @@ guchar GetVIS () {
readPcm(441);
// Apply Hann window
for (i = 0; i < 882; i++) in[i] = pcm.Buffer[pcm.WindowPtr + i - 441] / 32768.0 * Hann[i];
for (i = 0; i < 882; i++) fft.in[i] = pcm.Buffer[pcm.WindowPtr + i - 441] / 32768.0 * Hann[i];
// FFT of last 20 ms
fftw_execute(Plan2048);
@ -53,7 +53,7 @@ guchar GetVIS () {
// Find the bin with most power
MaxBin = 0;
for (i = 0; i <= GetBin(6000, FFTLen); i++) {
Power[i] = pow(out[i], 2) + pow(out[FFTLen - i], 2);
Power[i] = power(fft.out[i]);
if ( (i >= GetBin(500,FFTLen) && i < GetBin(3300,FFTLen)) &&
(MaxBin == 0 || Power[i] > Power[MaxBin]))
MaxBin = i;
@ -159,7 +159,6 @@ guchar GetVIS () {
if (++ptr == 25) {
setVU(Power, 2048, 6);
pcm.PeakVal = 0;
ptr = 0;
}