1
$\begingroup$

I have a problem with a FFT analysis in MATLAB, which is probably related to my limited understanding of the fundamentals of Fourier analysis.

  • I use MATLAB’s built-in fft analysis function.
  • The sampling frequency is 48000Hz The samples signal is a fairly clean sinus at 3000Hz (se below)
  • I use a flat top windowing.
  • I sample two voltage signals and use them to calculate the impedance of a device under test. An example of time data fed into the fft function, after windowing, is presented below enter image description here

When I increases the sampling time from shorter to a longer times, e.g. 0.5sec to 1sec, I get incorrect impedance results (e.g. negative resistance). Intuitively I would assume that longer sampling times would increase the accuracy. I have tracked the difference between the two different sample times to a change of the sign of the imaginary part of the discrete Fourier transform output for the 3000Hz-frequency bin. I do not know if this can be related to the incorrect results. For the shorter sampling time the imaginary part of the DFT is negative, se left figure below. And for the longer sampling time the imaginary part of the DFT is positive, se right figure below.

http://i.imgur.com/VOSPNgw.gif

I am grateful for any ideas of what might be causing this issue? I’m also happy to provide additional information if needed.

Kind regards, Johan

edit: My code is quite long, but below I have attached a code extract which I think is most relevant to the problem. I have not written all this code myself, but the author is unfortunately not available to me for questions. Generally this code have worked very well for a long time, the issue described in my original post is the only time that I have ever seen a problem.

% stSet.freq is the excitation frequency
% TH is the time data matrix

[Gij, freq, spec] = specDens(TH(inx,2:end), measData.Fs, blockSize, windowFcn, nOverlap, nFFT); 
inxHD = ([1:nHarmonics].*stSet.freq)./spec.x(2)+1; %% Defining indicies for HD components
inxHD = inxHD(inxHD<nFFT/2-1);%% Taking only HD components < Fs/2

% I = U2/R
% U = (U1-U2)
% Hiu = U/I = (U1-U2)/U2*R (electrical impedance)

% Gij      - cross- and autopower spectral density functions 

% Channels
ch_U1   = 1; % Ch.#1 - U1
ch_U2   = 2; % Ch.#2 - U2
R = 1.; % Resistance of measurement resistor

H.iu    = (Gij(inxHD,ch_U1,ch_U1) - Gij(inxHD,ch_U2,ch_U1)) ./ Gij(inxHD,ch_U2,ch_U1)*R; % Electrical impedance

Subfunction: specDens

function [Gij, freq, spec] = specDens(TH, Fs, blockSize, windowFcn, nOverlap, nFFT)
%
% Syntax:
% [Gij, freq, spec] = specDensRE2(TH, Fs, blockSize, windowFcn, nOverlap, nFFT)
%
% MATLAB function for computing cross- and autopower spectral density functions from continous signal measurements
%
%  __Outputs__
%   Gij      - cross- and autopower spectral density functions
%   freq     - frequency vector
%   spec     - Single sided linear spectrum
%
%  __Inputs__
%   TH        - measured signal, [m, n] = size(TH), m = No of samples, n = No of channels
%   Fs        - sampling frequency
%   blockSize - Block size
%   windowFcn - user defined window,  no input => Hanning
%   nOverlap  - user defined overlap, no input => nOverlap <=> 75%
%   nFFT      - Number of fft points


if(nargin<6); 
    nFFT = blockSize; end;

if(nargin<5); 
    nOverlap = round(blockSize*.75); end;

if(nargin<4); 
    [windowFcn, windowLSS] = window_re(blockSize,'Hanning'); end;

if(ischar(windowFcn)); %% assuming windowFcn is defining window type
    [windowFcn, windowLSS] = window_re(blockSize,windowFcn); end;


[nSamples, nCh] = size(TH); %% nSamples=No of samples, nCh=No of channels
if(nSamples<nCh);
    disp('Assuming input data transposed.');
    TH = TH.';
    [nSamples, nCh] = size(TH); %% nSamples=No of samples, nCh=No of channels
end

nAvg        = fix((nSamples-nOverlap)/(blockSize-nOverlap));
nLost       = mod((nSamples-nOverlap),(blockSize-nOverlap));
windowFcn   = windowFcn(:);                             %% making sure windowFcn is a column vector
Gij         = zeros(nFFT,nCh,nCh);                      %% cross- and autopower spectral density functions
scaleFactor = (nFFT/blockSize)*nAvg*sum(windowFcn)^2;   %% quad. normalising scale factor
freq        = [0:(nFFT-1)]'*Fs/nFFT;                    %% frequency vector
tmp         = zeros(nFFT,nCh);

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% % % CROSS- AND POWER SPECTRA
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
inxSelect = [1:blockSize];
%disp(['Calculating spectral density functions, nAvg = ' int2str(nAvg) ', nLost = ' int2str(nLost)]);
for loop=1:nAvg
    tmp = TH(inxSelect,:)-repmat(mean(TH(inxSelect,:)),blockSize,1); %% Removing BIAS
    %tmp = detrend(TH(inxSelect,:));     %% select a block and de-trend

    tmp = tmp.*repmat(windowFcn,1,nCh);  %% apply window to signal(s)
    tmp = fft(tmp,nFFT);                 %% FFT 
    for ii = 1:nCh
        for jj = ii:nCh
            Gij(:,ii,jj) = Gij(:,ii,jj) + tmp(:,ii).*conj(tmp(:,jj));
        end
    end
    inxSelect = inxSelect + round(blockSize-nOverlap);
end
for ii = 1:nCh
    for jj = (ii+1):nCh
        Gij(:,jj,ii) = conj(Gij(:,ii,jj));
    end
end
Gij = Gij/scaleFactor;


% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% % % IF NUMBER OF OUTPUT ARGUMENTS > 2
% % %    => OUTPUT SINGLE SIDED AMPLITUDE SPECTRUM
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
if(nargout>2);
    % % Selecting single sided spectrum (assuming real signal)
    if(rem(nFFT,2)==1);    % nFFT odd
        inx = [1:(nFFT+1)/2]';
    else
        inx = [1:nFFT/2+1]';
    end
    spec.x  = freq(inx);       %% frequency vector
    spec.y  = zeros(length(inx),nCh);
    for ii=1 : nCh
        spec.y(:,ii)=sqrt(Gij(inx,ii,ii)*2);
    end

    %% PSD
    % PSD       = (Gij(inx,ii,ii)*2)/df/windowLSS; %% needs window LSS to calc from AP 
    % PSD       = spec.y .^2/df/windowLSS;        
    % RMS_PSD   = sqrt(df*sum(PSD));  

    spec.Fs        = Fs;    
    spec.nFFT      = nFFT;
    spec.RMSinfo   = char('RMS(time) = sqrt(mean(abs(x).^2))','RMS(spec) = sqrt(sum(spec.y.^2)/windowLSS)','NOK for signals with BIAS ~= 0');
end
$\endgroup$
-1
$\begingroup$

actually decreasing sampling frequency may create the aliasing effect (sampled signal will be distorted and nothing similar to the original one)

http://en.wikipedia.org/wiki/Aliasing

if You want to properly sample Your signal use Nyquist - Shannon theory, which states that sampling frequency should be at least twice the signal frequency (it was proven that twice is enough and further increasing sampling rate does not gives better results)

http://en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem

thus Your sampling rate should be at least 6000Hz (sampling frequency of the sinus signal)

could You also provide the matlab code You used in this?

$\endgroup$
  • $\begingroup$ Thanks for the answer. The same sampling frequency 48kHz is used in both cases. My signal frequency is 3kHz; this gives 16 samples per period which I think should be sufficient. I also analyze the harmonic distortion, which requires me to use a higher sampling frequency than the Nyquist frequency. Please see the code which I added to the original post. $\endgroup$ – JohanGustafsson Feb 15 '14 at 21:52
  • $\begingroup$ to be honest I was confused when You mentioned the "When I increases the sampling time from shorter to a longer times, e.g. 0.5sec to 1sec" the code for the fft itself looks for me ok $\endgroup$ – akfaz Feb 15 '14 at 22:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.