MATLAB: To Get Correct FFT Result the length(x) must be An Integral Multiple of 100

fftMATLAB

Hi, all

I'm trying to get the amplitude- frequency relationship and phase- frenquency relationship using fft. My signal is given by:

t=0:0.005:14; % Time [s] Where mod(length(t),100)~=0
% t=t=0:0.005:14; t(end)=[]; % Where mod(length(t),100)==0
Vs=10*cos(2*pi*20*t-pi/5)+13*cos(2*pi*60*t+pi/6);
Spectrum=fft(Vs)/length(Vs);
Spectrum(2:end-1)=Spectrum(2:end-1)*2;

To get the amplitude and the phase of the signal of specified frequency, I used:

Amplitude=abs(Spectrum);
Phase=angle(Spectrum);

The question is, if the length(Vs) is an integral multiple of 100 (100,2700, etc.), the programe returns a correct amplitude and phase like:

Both Amplitude and Phase value at the DC, 20Hz, 40Hz is the setting value. However, once the length(Vs) is not of an integral multiple of 100, the same program returns a totally wrong result like:

I don't know the reason, and I want a method by which I can get right Spectrum no matter what the length(Vs) is.(On the premise of guaranteeing sampling theorem).

Thanks

Best Answer

The desired plot will be obtained when N, the number of elements of t, is an integer multiple of the ratios Fs/Fn, where Fs is the sample frequency and Fn are the frequencies in Vs (in Hz). In this case we have

Fs = 200;
Fn = [20 60];

Given these vaues, selecting N = k*100 will satisy the criterion:

syms k N
N = k*sym(100);
N*sym(Fn)./sym(Fs)  % integer valued for all integer k
ans =

Here is the "good" case in the question:

N = 2800;
N.*Fn./Fs   % integer

ans = 1×2
   280   840
tfinal = func(N,Fs,Fn)
tfinal = 13.9950

Here is another case, where N satisfies the criteria, but is not a muliple of 100

N = 2820;
N.*Fn./Fs   % integer
ans = 1×2
   282   846
tfinal = func(N,Fs,Fn)
tfinal = 14.0950

As you already showed, choosing N = 2801, which does not satisfy the criterion, results in an undesired plot

N = 2801;
N.*Fn./Fs  % not both integers
ans = 1×2
  280.1000  840.3000
tfinal = func(N,Fs,Fn)
tfinal = 14

Note this result is not wrong. It is correct based on the input parameters.

function tfinal = func(N,Fs,Fn)
t = (0:(N-1))/Fs;
Vs = 10*cos(2*pi*Fn(1)*t-pi/5)+13*cos(2*pi*Fn(2)*t+pi/6);
Spectrum = fft(Vs)/length(Vs);
Spectrum(2:end-1) = Spectrum(2:end-1)*2;
Amplitude = abs(Spectrum);
Phase = angle(Spectrum);
figure;
stem((0:N-1)/N*Fs,Amplitude)
xlim([0 Fs/2])
xlabel('Hz')
tfinal = t(end);
end

Related Solutions

MATLAB: I get the peaks at the wrong frequency and with the wrong amplitude.

hello Anna

look at my fft function at the end of my code

clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% data
data = importdata('SFC5S_nov25_ST_1_1.txt');
dt = 1e-6; % 1 micro seconds
signal = data(:,1);
samples = length(signal);
Fs = 1/dt; % sampling frequency (Hz)
%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 1;
if decim>1
    signal = decimate(signal,decim);
    Fs = Fs/decim;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 512;    % 
OVERLAP = 0.95;
% spectrogram dB scale
spectrogram_dB_scale = 80;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% options 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if you are dealing with acoustics, you may wish to have A weighted
% spectrums 
% option_w = 0 : linear spectrum (no weighting dB (L) )
% option_w = 1 : A weighted spectrum (dB (A) )
option_w = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% display 1 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[freq, sensor_spectrum] = myfft_peak(signal,Fs,NFFT,OVERLAP);
% convert to dB scale (ref = 1)
sensor_spectrum_dB = 20*log10(sensor_spectrum);
% apply A weigthing if needed

if option_w == 1
    pondA_dB = pondA_function(freq);
    sensor_spectrum_dB = sensor_spectrum_dB+pondA_dB;
    my_ylabel = ('Amplitude (dB (A))');
else
    my_ylabel = ('Amplitude (dB (L))');
end
figure(1),plot(freq,sensor_spectrum_dB,'b');grid
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
xlabel('Frequency (Hz)');ylabel(my_ylabel);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% display 2 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[sg,fsg,tsg] = specgram(signal,NFFT,Fs,hanning(NFFT),floor(NFFT*OVERLAP));  
% FFT normalisation and conversion amplitude from linear to dB (peak)
sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg));     % NB : X=fft(x.*hanning(N))*4/N; % hanning only
% apply A weigthing if needed
if option_w == 1
    pondA_dB = pondA_function(fsg);
    sg_dBpeak = sg_dBpeak+(pondA_dB*ones(1,size(sg_dBpeak,2)));
    my_title = ('Spectrogram (dB (A))');
else
    my_title = ('Spectrogram (dB (L))');
end
% saturation of the dB range : 
% saturation_dB = 60;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;
sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;
% plots spectrogram
figure(2);
imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
axis('xy');colorbar('vert');grid
title([my_title ' / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Frequency (Hz)');
function pondA_dB = pondA_function(f)
	% dB (A) weighting curve
	n = ((12200^2*f.^4)./((f.^2+20.6^2).*(f.^2+12200^2).*sqrt(f.^2+107.7^2).*sqrt(f.^2+737.9^2)));
	r = ((12200^2*1000.^4)./((1000.^2+20.6^2).*(1000.^2+12200^2).*sqrt(1000.^2+107.7^2).*sqrt(1000.^2+737.9^2))) * ones(size(f));
	pondA = n./r;
	pondA_dB = 20*log10(pondA(:));
end
function  [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal  (example sinus amplitude 1   = 0 dB after fft).
% Linear averaging
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer overlap % (between 0 and 0.95)
signal = signal(:);
samples = length(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,1);
    s_tmp((1:samples)) = signal;
    signal = s_tmp;
    samples = nfft;
end
% window : hanning
window = hanning(nfft);
window = window(:);
%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows
    % main loop
    fft_spectrum = 0;
    for i=1:spectnum
        start = (i-1)*offset;
        sw = signal((1+start):(start+nfft)).*window;
        fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft);     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end
    fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling
% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_spectrum = fft_spectrum(select);
freq_vector = (select - 1)*Fs/nfft;
end

MATLAB: How can i find the peaks of the signal and get a phase and magnitude plot

The easiest way is to create and plot a one-sided fft (link). The phase plot uses the angle function.

To find the peaks and their locations, use the Signal Processing Toolbox findpeaks (link) function.

Example —

N = 1000;
t = linspace(0, 10*pi, N);                      % Time Vector
Fs = 1/(t(2)-t(1));                             % Sampling Frequency
Fn = Fs/2;                                      % Nyquist Frequency
s = sin(2*pi*5*t);                              % Signal Vector
FTs = fft(s)/N;                                 % Divide By ‘N’ To Normalise Amplitudes
Fv = linspace(0, 1, fix(N/2)+1)*Fn;             % Frequency Vector (1-Sided)
Iv = 1:length(Fv);                              % Index Vector
figure(1)
subplot(2,1,1)
plot(Fv, abs(FTs(Iv))*2)
grid
ylabel('Amplitude')
title('Amplitude Plot')
subplot(2,1,2)
plot(Fv, angle(FTs(Iv))*180/pi)
grid
xlabel('Frequency (Hz)')
ylabel('Phase (°)')
title('Phase Plot')

Best Answer

Related Solutions

MATLAB: I get the peaks at the wrong frequency and with the wrong amplitude.

MATLAB: How can i find the peaks of the signal and get a phase and magnitude plot

Related Question