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
