# FFT analysis of longer data sets

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 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. 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


## migrated from stackoverflow.comMar 21 '14 at 14:19

This question came from our site for professional and enthusiast programmers.

## 1 Answer

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?

• 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. – JohanGustafsson Feb 15 '14 at 21:52
• 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 – akfaz Feb 15 '14 at 22:01