Trying to debug my FFT

[VelekHacks]

I have my modified version of the RX app trying to sample from the 2m Ham band, and I’m having trouble getting things to match up.
My tinySA Spectrum Analyzer sees the right peak.

/**
    @file   singleRX.cpp
    @author Lime Microsystems (www.limemicro.com)
    @brief  RX example
 */
// Try to get this to demodulate FM radio with the blog entry at:
// https://witestlab.poly.edu/blog/capture-and-decode-fm-radio/
// #define USE_GNU_PLOT
#include "lime/LimeSuite.h"
#include <iostream>
#include <chrono>
#include <fstream>
#ifdef USE_GNU_PLOT
#include "gnuPlotPipe.h"
#endif
#include <complex>
#include <thread>
#include <queue>
#include <mutex>
#include <condition_variable>
#include <boost/lockfree/queue.hpp>
#include <fftw3.h>
using namespace std;

lms_device_t *device = NULL;
std::mutex mtx;
boost::lockfree::queue<float *> sampleQueue(128);
std::atomic<bool> dataAvailable(false);
std::atomic<bool> exitFlag(false);
std::condition_variable cv;

const double bandStart = 144e6;                                      // Start from the Nyquist frequency
const double bandEnd = 148e6;                                        // End at the Nyquist frequency
const double startFrequency = bandStart + (bandEnd - bandStart) / 2; // Start at the center of the FM band
const int bandWidth = bandEnd - bandStart;
const int sampleRate = bandWidth * 2;
const int bufferSize = 4 * 1024; // complex samples per buffer

std::vector<double> calculatePower(const std::vector<float> &samples)
{
    // Calculate the number of channels
    int binCount = samples.size() / 2;

    // Create a vector to hold the power of each channel
    std::vector<double> power(binCount);

    // Create input array
    fftw_complex *input = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * samples.size() / 2);

    int n = 0;
    // Copy the samples for this channel to the input array
    for (long unsigned int i = 0; i < samples.size(); i += 2)
    {
        input[n][0] = samples[i];
        input[n][1] = samples[i + 1];
        n++;
    }
    auto start = std::chrono::high_resolution_clock::now();

    // Create the FFTW plan
    fftw_plan plan = fftw_plan_dft_1d(binCount, input, input, FFTW_FORWARD, FFTW_ESTIMATE);

    // Execute the FFT
    fftw_execute(plan);
    auto end = std::chrono::high_resolution_clock::now();

    // Calculate the power of the signal for this channel
    for (int i = 0; i < binCount; i++)
    {
        double real = input[i][0];
        double imag = input[i][1];
        power[i] += std::round(real * real + imag * imag);
    }

    // Clean up
    fftw_destroy_plan(plan);
    fftw_free(input);

    // Calculate and print the elapsed time in microseconds
    auto duration = std::chrono::duration_cast<std::chrono::microseconds>(end - start);
    std::cout << "FFT took " << duration.count() << " microseconds\n";

    for (long unsigned int i = 0; i < power.size(); i++)
    {
        double freq = startFrequency + i * bandWidth / power.size(); // Subtract the width of each band
        if (power[i] > 100000)
        {
            cout << "Band " << i << " starts at " << freq / 1e6 << " MHz and has power " << power[i] << "\n";
        }
    }

    return power;
}

int error()
{
    if (device != NULL)
        LMS_Close(device);
    exit(-1);
}

void waterfallThread()
{
    while (true)
    {
        float *buffer = nullptr;
        while (!sampleQueue.pop(buffer))
        {
            // Wait for data to become available
            std::unique_lock<std::mutex> lock(mtx);
            cv.wait(lock, []
                    { return dataAvailable.load() || exitFlag.load(); }); // Check the flag here
            dataAvailable.store(false);
        }
        cout << "Got buffer\n";

        std::vector<double> power = calculatePower(std::vector<float>(buffer, buffer + bufferSize));

        double binWidth = bandWidth / power.size();
        cout << "Bin width: " << binWidth << "\n";
        for (long unsigned int i = 0; i < power.size(); i++)
        {
            double freq = bandStart + i * binWidth; // Subtract the width of each band
            if (power[i] > 1000)
            {
                cout << "Band " << i << " starts at " << freq / 1e6 << " MHz and has power " << power[i] << "\n";
            }
        }

        if (exitFlag.load())
            break;

        // Don't forget to delete the buffer when you're done with it
        delete[] buffer;
    }
}

int main(int argc, char **argv)
{
    // Find devices
    // First we find number of devices, then allocate large enough list,  and then populate the list
    int n;

    if ((n = LMS_GetDeviceList(NULL)) < 0) // Pass NULL to only obtain number of devices
        error();
    cout << "Devices found: " << n << endl;
    if (n < 1)
        return -1;

    lms_info_str_t *list = new lms_info_str_t[n]; // allocate device list

    if (LMS_GetDeviceList(list) < 0) // Populate device list
        error();

    for (int i = 0; i < n; i++) // print device list
        cout << i << ": " << list[i] << endl;
    cout << endl;

    // Open the first device
    if (LMS_Open(&device, list[0], NULL))
        error();

    delete[] list; // free device list

    // Initialize device with default configuration
    // Do not use if you want to keep existing configuration
    // Use LMS_LoadConfig(device, "/path/to/file.ini") to load config from INI
    if (LMS_Init(device) != 0)
        error();

    // Enable RX channel
    // Channels are numbered starting at 0
    if (LMS_EnableChannel(device, LMS_CH_RX, 0, true) != 0)
        error();

    if (LMS_EnableChannel(device, LMS_CH_TX, 0, true) != 0) // Fix for v2
        error();

    // Set center frequency to 800 MHz
    // int F_offset = 250000;
    // int desired_freq = 98 * 1e5 - F_offset;
    int desired_freq = startFrequency;
    if (LMS_SetLOFrequency(device, LMS_CH_RX, 0, desired_freq) != 0)
        error();

    // print currently set center frequency
    float_type freq;
    if (LMS_GetLOFrequency(device, LMS_CH_RX, 0, &freq) != 0)
        error();
    cout << "\nCenter frequency: " << freq / 1e6 << " MHz\n";

    // select antenna port
    lms_name_t antenna_list[10]; // large enough list for antenna names.
                                 // Alternatively, NULL can be passed to LMS_GetAntennaList() to obtain number of antennae
    if ((n = LMS_GetAntennaList(device, LMS_CH_RX, 0, antenna_list)) < 0)
        error();

    cout << "Available antennae:\n"; // print available antennae names
    for (int i = 0; i < n; i++)
        cout << i << ": " << antenna_list[i] << endl;

    if ((n = LMS_GetAntenna(device, LMS_CH_RX, 0)) < 0) // get currently selected antenna index
        error();
    // print antenna index and name
    cout << "Automatically selected antenna: " << n << ": " << antenna_list[n] << endl;

    if (LMS_SetAntenna(device, LMS_CH_RX, 0, LMS_PATH_LNAH) != 0) // manually select antenna
        error();

    if ((n = LMS_GetAntenna(device, LMS_CH_RX, 0)) < 0) // get currently selected antenna index
        error();
    // print antenna index and name
    cout << "Manually selected antenna: " << n << ": " << antenna_list[n] << endl;

    // Set sample rate to 8 MHz, preferred oversampling in RF 8x
    // This set sampling rate for all channels
    if (LMS_SetSampleRate(device, sampleRate, 8) != 0)
        error();
    // print resulting sampling rates (interface to host , and ADC)
    float_type rate, rf_rate;
    if (LMS_GetSampleRate(device, LMS_CH_RX, 0, &rate, &rf_rate) != 0) // NULL can be passed
        error();
    cout << "\nHost interface sample rate: " << rate / 1e6 << " MHz\nRF ADC sample rate: " << rf_rate / 1e6 << "MHz\n\n";

    // Example of getting allowed parameter value range
    // There are also functions to get other parameter ranges (check LimeSuite.h)

    // Get allowed LPF bandwidth range
    lms_range_t range;
    if (LMS_GetLPFBWRange(device, LMS_CH_RX, &range) != 0)
        error();

    cout << "RX LPF bandwidth range: " << range.min / 1e6 << " - " << range.max / 1e6 << " MHz\n\n";

    // Configure LPF, bandwidth 8 MHz
    if (LMS_SetLPFBW(device, LMS_CH_RX, 0, sampleRate / 2) != 0)
        error();

    // Set RX gain
    if (LMS_SetNormalizedGain(device, LMS_CH_RX, 0, 0.7) != 0)
        error();
    // Print RX gain
    float_type gain; // normalized gain
    if (LMS_GetNormalizedGain(device, LMS_CH_RX, 0, &gain) != 0)
        error();
    cout << "Normalized RX Gain: " << gain << endl;

    unsigned int gaindB; // gain in dB
    if (LMS_GetGaindB(device, LMS_CH_RX, 0, &gaindB) != 0)
        error();
    cout << "RX Gain: " << gaindB << " dB" << endl;

    // Perform automatic calibration
    if (LMS_Calibrate(device, LMS_CH_RX, 0, sampleRate, 0) != 0)
        error();

    // Enable test signal generation
    // To receive data from RF, remove this line or change signal to LMS_TESTSIG_NONE
    //    if (LMS_SetTestSignal(device, LMS_CH_RX, 0, LMS_TESTSIG_NCODIV8, 0, 0) != 0)
    //        error();

    // Streaming Setup

    // Initialize stream
    lms_stream_t streamId;
    streamId.channel = 0;                         // channel number
    streamId.fifoSize = bufferSize * bufferSize;  // fifo size in samples
    streamId.throughputVsLatency = 1.0;           // optimize for max throughput
    streamId.isTx = false;                        // RX channel
    streamId.dataFmt = lms_stream_t::LMS_FMT_F32; // 32-bit floats
    if (LMS_SetupStream(device, &streamId) != 0)
        error();

    LMS_StartStream(&streamId);

#ifdef USE_GNU_PLOT
    GNUPlotPipe gp;
    gp.write("set size square\n set xrange[-1:1]\n set yrange[-1:1]\n");
#endif
    auto t1 = chrono::high_resolution_clock::now();
    auto t2 = t1;

    // In your main function or wherever you're receiving the samples:
    string filename = "samples_" + to_string(freq / 1e6) + ".csv";
    // std::thread writerThread(writeSamplesToFile, filename);
    std::thread writerThread(waterfallThread);

    // long long timeForBuffer = (1 / sampleRate) * bufferSize * 1e6; // time in microseconds

    while (chrono::high_resolution_clock::now() - t1 < chrono::milliseconds(5)) // run for 10 seconds
    {
        int samplesRead;
        // Allocate a new buffer
        float *buffer = new float[bufferSize * 2];

        lms_stream_status_t status;
        if (LMS_GetStreamStatus(&streamId, &status) != 0)
        {
            // Handle error
            cout << "Stream Status error" << endl;
        }

        // Receive samples
        samplesRead = LMS_RecvStream(&streamId, buffer, bufferSize, NULL, 1000);
        if (samplesRead <= 0)
        {
            cout << "error: " << samplesRead << endl;
        }
        cout << "Read " << samplesRead << " samples.\n";
        while (!sampleQueue.push(buffer))
            ;

        dataAvailable.store(true);
        cv.notify_one();

        // I and Q samples are interleaved in buffer: IQIQIQ...
        /*
            INSERT CODE FOR PROCESSING RECEIVED SAMPLES
        */

        // Print stats (once per second)
        if (chrono::high_resolution_clock::now() - t2 > chrono::seconds(1))
        {
            t2 = chrono::high_resolution_clock::now();
            lms_stream_status_t status;
            // Get stream status
            LMS_GetStreamStatus(&streamId, &status);
            cout << "RX data rate: " << status.linkRate / 1e6 << " MB/s\n";                       // link data rate
            cout << "RX fifo: " << 100 * status.fifoFilledCount / status.fifoSize << "%" << endl; // percentage of FIFO filled
        }
    }

    cout << "Finished sampling" << endl;

    // Stop streaming
    LMS_StopStream(&streamId);            // stream is stopped but can be started again with LMS_StartStream()
    LMS_DestroyStream(device, &streamId); // stream is deallocated and can no longer be used

    // Close device
    LMS_Close(device);

    cout << "Waiting for samples to be written to file" << endl;

    while (!sampleQueue.empty())
    {
        std::this_thread::sleep_for(std::chrono::milliseconds(100));
    }

    cout << "Complete" << endl;

    // Don't forget to join the thread at the end of your program:
    exitFlag.store(true);
    cv.notify_one();
    writerThread.join();

    return 0;
}

When I key the mike on my Yaesu on 145MHz, I get this output from the program:

Band 0 starts at 144 MHz and has power 1077
Band 107 starts at 144.209 MHz and has power 1372
Band 112 starts at 144.219 MHz and has power 2117
Band 129 starts at 144.252 MHz and has power 1121
Band 130 starts at 144.254 MHz and has power 2473
Band 135 starts at 144.264 MHz and has power 1448
Band 136 starts at 144.266 MHz and has power 1108
Band 137 starts at 144.268 MHz and has power 1137
Band 138 starts at 144.27 MHz and has power 1203
Band 139 starts at 144.271 MHz and has power 1205
Band 140 starts at 144.273 MHz and has power 5078
Band 141 starts at 144.275 MHz and has power 2566
Band 142 starts at 144.277 MHz and has power 2350
Band 143 starts at 144.279 MHz and has power 1734
Band 145 starts at 144.283 MHz and has power 16861
Band 146 starts at 144.285 MHz and has power 12249
Band 147 starts at 144.287 MHz and has power 13518
Band 148 starts at 144.289 MHz and has power 17202
Band 149 starts at 144.291 MHz and has power 8001
Band 150 starts at 144.293 MHz and has power 38635
Band 151 starts at 144.295 MHz and has power 72927
Band 152 starts at 144.297 MHz and has power 187574
Band 153 starts at 144.299 MHz and has power 1.23797e+06
Band 154 starts at 144.301 MHz and has power 3.54064e+06
Band 155 starts at 144.303 MHz and has power 260385
Band 156 starts at 144.305 MHz and has power 90073
Band 157 starts at 144.307 MHz and has power 48924
Band 158 starts at 144.309 MHz and has power 69522
Band 159 starts at 144.311 MHz and has power 9174
Band 160 starts at 144.312 MHz and has power 8282
Band 161 starts at 144.314 MHz and has power 6145
Band 162 starts at 144.316 MHz and has power 3649
Band 163 starts at 144.318 MHz and has power 31159
Band 164 starts at 144.32 MHz and has power 4938
Band 165 starts at 144.322 MHz and has power 3310
Band 166 starts at 144.324 MHz and has power 2349
Band 167 starts at 144.326 MHz and has power 1152
Band 168 starts at 144.328 MHz and has power 6467
Band 169 starts at 144.33 MHz and has power 2921
Band 170 starts at 144.332 MHz and has power 2108
Band 171 starts at 144.334 MHz and has power 1687
Band 172 starts at 144.336 MHz and has power 1405
Band 173 starts at 144.338 MHz and has power 2480
Band 174 starts at 144.34 MHz and has power 1700
Band 175 starts at 144.342 MHz and has power 1461
Band 176 starts at 144.344 MHz and has power 1831
Band 178 starts at 144.348 MHz and has power 1040
Band 179 starts at 144.35 MHz and has power 1084
Band 180 starts at 144.352 MHz and has power 1190
Band 181 starts at 144.353 MHz and has power 3003
Band 186 starts at 144.363 MHz and has power 1150
Band 191 starts at 144.373 MHz and has power 1027
Band 196 starts at 144.383 MHz and has power 1058
Band 214 starts at 144.418 MHz and has power 1021
Band 759 starts at 145.482 MHz and has power 2543
Band 763 starts at 145.49 MHz and has power 4453
Band 764 starts at 145.492 MHz and has power 4017
Band 768 starts at 145.5 MHz and has power 10502
Band 772 starts at 145.508 MHz and has power 1019
Band 773 starts at 145.51 MHz and has power 7152
Band 777 starts at 145.517 MHz and has power 1639
Band 778 starts at 145.519 MHz and has power 1212
Band 973 starts at 145.9 MHz and has power 1246
Band 977 starts at 145.908 MHz and has power 1029
Band 1564 starts at 147.054 MHz and has power 1523
Band 1573 starts at 147.072 MHz and has power 1120
Band 1578 starts at 147.082 MHz and has power 6582
Band 1582 starts at 147.09 MHz and has power 7136
Band 1583 starts at 147.092 MHz and has power 7090
Band 1586 starts at 147.097 MHz and has power 2092
Band 1587 starts at 147.099 MHz and has power 194056
Band 1588 starts at 147.101 MHz and has power 4924
Band 1589 starts at 147.103 MHz and has power 1669
Band 1590 starts at 147.105 MHz and has power 1394
Band 1591 starts at 147.107 MHz and has power 2970
Band 1592 starts at 147.109 MHz and has power 10582
Band 1596 starts at 147.117 MHz and has power 3848
Band 1597 starts at 147.119 MHz and has power 1436

[RightHalfPlane]

Most curious when I run it I get -

Got buffer
FFT took 47 microseconds
Band 1792 starts at 149.5 MHz and has power 346433
Bin width: 1953
Band 1791 starts at 147.498 MHz and has power 1212
Band 1792 starts at 147.5 MHz and has power 346433
Band 1793 starts at 147.502 MHz and has power 1417
Waiting for samples to be written to file
Complete

[RightHalfPlane]

I took one of my simple SoapySDR programs and modified it to do the FFT calculations.

I widen the bandwidth to 6MHZ - to avoid the droop at the edges.

I also put in the terms before doing the FFT

    buf1[2*n] *= pow(-1.0,n);
    buf1[2*n+1] *= pow(-1.0,n);

to rotate the FFT to its desired orientation. My program returns the correct frequencies. Putting the rotation into your program it returns 145.5 instead of 145 - so there appears to be another problem with your program’s calculations - it should be closer than that. My program writes a I/Q file that you can open with SdrGlut or Gqrx to display the power spectrum. I have attached an example from a run. It shows the expected 145 MHZ spike plus another spike at 146 MHZ that the LimeMini 2.0 is transmitting for no reason.

#include <SoapySDR/Version.hpp>
#include <SoapySDR/Modules.hpp>
#include <SoapySDR/Registry.hpp>
#include <SoapySDR/Device.hpp>
#include <SoapySDR/Formats.h>
#include <stdio.h> //printf
#include <stdlib.h> //free
#include <complex.h>
#include <signal.h>
#include <fftw3.h>

// c++ -std=c++11 -o simpleFFT01 simpleFFT01.cpp -lfftw3f -Wall -lSoapySDR

int process(float *buf1,float *buf2,int length,fftwf_plan p1);

int loop = 1;
void sigIntHandler(const int)
{
    loop = 0;
}

double samplerate=6e6;

int main(void)
{
		
	signal(SIGINT, sigIntHandler);

    std::string argStr;
        
    std::vector<SoapySDR::Kwargs> results;
    
    results = SoapySDR::Device::enumerate();
    
    printf("Number of Devices Found: %ld\n",(long)results.size());
    
    if(results.size() < 1)return 1;

    SoapySDR::Kwargs deviceArgs;

   for(unsigned int k=0;k<results.size();++k){
    		printf("SDR device =  %ld ",(long)k);
			deviceArgs = results[k];
			for (SoapySDR::Kwargs::const_iterator it = deviceArgs.begin(); it != deviceArgs.end(); ++it) {
				if (it->first == "label")printf(" %s = %s\n ",it->first.c_str(), it->second.c_str());
			}
    }
    
    printf("\n");

    SoapySDR::Device *sdr = SoapySDR::Device::make(results[0]);

	if (sdr == NULL)
	{
		printf("SoapySDRDevice_make fail: \n");
		return EXIT_FAILURE;
	}
	
	sdr->setSampleRate(SOAPY_SDR_RX, 0, samplerate);
			
	sdr->setFrequency(SOAPY_SDR_RX, 0, 146e6);
	
			
	const std::vector<size_t> channels = {(size_t)0};
	
	SoapySDR::Stream *rxStream;
	
	try {
		rxStream = sdr->setupStream(SOAPY_SDR_RX, SOAPY_SDR_CF32, channels);
    } catch(const std::exception &e) {
        std::string streamExceptionStr = e.what();
        printf("Error: %s\n",streamExceptionStr.c_str());
        exit(1);
    }

	sdr->activateStream(rxStream, 0, 0, 0); 
	
	sdr->setGain(SOAPY_SDR_RX, 0, -12.0); 
	
	fprintf(stderr, "start\n");
	
	int length=8192;

	float buf1[8192*2],buf2[8192*2];
	
	fftwf_plan p1;
				
	p1 = fftwf_plan_dft_1d(length,(fftwf_complex *)buf1, (fftwf_complex *)buf2, FFTW_FORWARD, FFTW_ESTIMATE);
	
	FILE *out11;
	out11=fopen("test1_IQ_146000000_6000000_fc.raw","wb");
	if(!out11)exit(1);
 
	//receive some samples
	while (loop)
	{
		void *buffs[] = {buf1}; //array of buffers
		float *out;
		int flags=0; //flags set by receive operation
		long long timeNs=0; //timestamp for receive buffer
		int rec=length;
		out=buf1;
		while(rec > 0 && loop){
		    buffs[0]=out;
			int ret = sdr->readStream(rxStream, buffs, rec, flags, timeNs, 1000000L);
			//printf("ret=%d, flags=%d, timeNs=%lld rec %d out %f count %lld\n", ret, flags, timeNs, rec,out[0],countPrint++);
			if(ret < 0)fprintf(stderr,"Read Error\n");
			if(ret > 0){
			    fwrite(out,8,ret,out11);
				out += 2*(long long)ret;
				rec -= ret;
			}
		}
		process(buf1,buf2,length,p1);
	}
	
	sdr->deactivateStream(rxStream, 0, 0);

	sdr->closeStream(rxStream);

	SoapySDR::Device::unmake(sdr);

	if(out11)fclose(out11);

	printf("Done\n");
	return EXIT_SUCCESS;
}

int process(float *buf1,float *buf2,int length,fftwf_plan p1)
{

    for(int n=0;n<length;++n){
        buf1[2*n] *= pow(-1.0,n);
        buf1[2*n+1] *= pow(-1.0,n);
    }

	fftwf_execute(p1);
	
	float amax=0;
	int target=0;
	for(int n=0;n<length;++n){
		double v=buf2[2*n]*buf2[2*n]+buf2[2*n+1]*buf2[2*n+1];
		if(v > 0.0)v=sqrt(v);
		if(v > amax){
			amax=v;
			target=n;
		}
	}
    fprintf(stderr,"amax %g target %d frequency %g\n",amax,target,(146.0e6-0.5*samplerate)+samplerate*target/(double)length);
	return 0;
}