CODEWITHSUNDEEP
matlabsignal-processingfft
Paul
Paul

Reputation: 307

Shifting the FFT of a signal in Matlab

I'm trying to implement a spectral correlation function to plot with the surf function. I think I understand the idea of the SCF as described in a paper I read, but I'm having trouble implementing my function in Matlab. I've been following these instructions:

Pseudo-code

I'm mostly having trouble shifting my pieces of the data properly. Is there an easy way to achieve step 3?

Here's what I tried in my code: 

function [output] = spectral(x, N)
    % This function does cyclostationary spectral analysis 
    % on a data set and returns some features

    t = length(x); 
    samplesPerFrame = floor(t / N); 

    count = 1; 

    for alpha = -1:0.01:1



        % Split up the samples into frames
        % Have to leave some samples out if unevenly split
        for i = 1:N+1
            frange = ((i - 1) * samplesPerFrame + 1):(i * samplesPerFrame);
            if i == N+1
                break;
            end
            xFrame(i, :) = x(frange);
            ts = [1:length(xFrame(i,:))]; 
            shiftLeft = fft(xFrame(i, :) .* exp(-1 * 2 * pi * 1i * (alpha / 2) .* ts));
            shiftRight = fft(xFrame(i, :).* exp(2 * pi * 1i * (alpha / 2) .* ts)); 

            S(i,:) = (1 / samplesPerFrame) .* shiftLeft .* conj(shiftRight);

        end

        Savg(count, :) = mean(S, 1);
        Ssmooth(count, :) = smooth(Savg(count,:), 'moving'); 
        count = count + 1; 
    end
    output = Ssmooth; 
end

Upvotes: 0

Views: 1520

Answers (2)

Chad Spooner
Chad Spooner

Reputation: 39

The method of spectral correlation estimation you are attempting is something I refer to as the Time-Smoothing Method of spectral correlation estimation, or the TSM. The code you posted cannot provide the correct answer except in some trivial cases such as alpha = 0. The reason is that you need to adjust the cyclic periodogram for each frame by a complex phase factor to compensate for the fact that each data block is a delayed version of the one preceding it.

If you replace the line

S(i,:) = (1 / samplesPerFrame) .* shiftLeft .* conj(shiftRight);

with the two lines

S(i,:) = (1 / samplesPerFrame) .* shiftLeft .* conj(shiftRight);

S(i, :) = S(i, :) * exp(-1i * 2 * pi * alpha * i * samplesPerFrame);

you'll be able to estimate the SCF. I confirmed this by applying your original code and the modified code to a BPSK signal with bit rate (normalized) of 1/10. In this case, one of your alpha values in the loop over alpha will exactly coincide with the true cycle frequency of 1/10. Only the modified code gives the correct SCF for the bit-rate cycle frequency.

Please see my blog cyclostationary.wordpress.com for more detail and examples. In particular, I have a post on the TSM at http://cyclostationary.blog/2015/12/18/csp-estimators-the-time-smoothing-method. (Corrected this link 5/2/17.)

Upvotes: 1

lennon310
lennon310

Reputation: 12689

It looks good actually.

You may also try circshift(fft(xFrame(i, :)),[1,a]) to achieve shiftRight, and circshift(fft(xFrame(i, :)),[1,-a]) to get shiftLeft. Please note here a is integer, indicates the elements in xFrame(i, :) that you wish to move, and corresponds to Fs*a in frequency domain where Fs is your sampling rate.

Upvotes: 1

Related Questions