Estimation of impedance in MATLAB
The FFT is simply an efficient method for computing the DFT (Discrete Fourrier Transform), which is a convolution between an N sample input waveform \$ x_n\$, and each of the N complex basis functions of the form \$ u_n=exp(i2\pi fn)\$, where f is the frequency and varies from 0 to N-1.

The niave, direct, brute force way of doing the DFT requires \$O(2N^2)\$ multiplies and adds. The FFT uses a clever factorisation to radically reduce the number of adds and multiplies required to \$O(2N \log_2(N))\$, worth doing if N is more than few dozen, and you need all the frequencies.

If you only need the convolution at a single frequency, or a small handful of frequencies, then doing it directly, doing just one frequency of the DFT, is quicker than using an FFT, which is an 'all or nothing' calculation.

For instance, to compute the complex gain v/i at a single frequency, you would start by effectively generating a complex signal at that frequency, and output the real component of it to the device under test. 'Effectively' means that you do not actually need to explicitly generate the full complex cosine+sine signal. If the period is a multiple of 4 samples, that means you can generate just the cosine waveform, and shift it by 90 degrees, or just index it shifted, to get the sine for free.

Once you have the received the \$v_n\$ waveform, you convolve it with the complex reference \$u_n\$. This has real part \$cos(2\pi f n)\$ and imaginary part \$i.sin(2\pi f n)\$. The convolution as done as $$realpart= \sum_{n=0}^{N-1}v_n.cos(2\pi fn)$$This is the part of v that's in phase with the stimulus signal. You get the imaginary, quadrature, part by doing the same summation having substituted sin() for cos() as $$imagpart= \sum_{n=0}^{N-1}v_n.sin(2\pi fn)$$You of course do not need to compute the cos(2pifn) when doing the sums, it reduces to indexing through your waveform one sample at a time!

In C, the algorithm is

float* v_ptr = v_buffer;  
float* real_ptr = &amp;cos_table[0];  
float* imag_ptr = &amp;cos_table[OFFSET_FOR_SINE_START];  
float real_sum=0, imag_sum=0;  
for (i=0; i&lt;N; i++){  
    real_sum += *real_ptr++ * *v_ptr;  
    imag_sum += *imag_ptr++ * *v_ptr++;
}


As a detail, the above assumes that cos_table is long enough so that sine can index it to the end after starting 1/4 cycle in. You would have to do something tricky if it was only N long.

As this computation is exactly the same, at that frequency, that the FFT does at all frequencies, it is subject to the same limitations of periodicity. That means if your waveform is not periodic in your chosen analysis vector length N, then your convolution results will be inaccurate. In the FFT this is manifested as 'spectral leakage'.

As in the FFT, there are two solutions to this. The first is to ensure that your N sample waveform always contains an integer number of cycles. This is easy when you are generating your own waveform! When the signal is not under your control, windowing is usually used prior to the FFT to reduce spectral leakage.