MATLAB: I wish to get the spectrum (Magnitude vs Frequency(Hz)) of a pure tone(1 kHz). However, the piece of code I’m employing to do so isn’t working. What am I doing wrong

fftfftshift

%fft of puretone
%the cell puretone{1,1} contains the puretone
  fft_puretone{1,1} = abs(fftshift(fft(puretone{1,1})));
plot(x_axis,fft_puretone{1,1})

However, I am getting this plot:

whereas I should get a peak at 1kHz. I think there is some problem in using fftshift() and fft() but I can't figure out what. There maybe some problem with setting up the axis which I've taken as

x_bins = [-N/2:N/2-1] %where N is the number of samples in the sound
x_axis = x_bins*fs/N;

Best Answer

Hello Komal,

The whole point of using fftshift is to put zero frequency at the center of the resulting frequency array. If you zero-pad fft_puretone on one side only, it messes up that symmetry. So fft_puretone needs to be zero-padded equally on each side.

An an N-point fft with time array spacing of delt and frequency array spacing of delf always obeys delt*delf = 1/N. That remains true after subsequent zero-padding.

For the correct frequency domain amplitude, you need to multiply the fft by a factor of 1/N.

You did not specify what puretone looks like in the time domain, but I am speculating that the time record was 5 seconds, and the amplitude of the wave = 1/sqrt(2) (peak), = 1/2 (rms).

Here is an example of how this would work. The resulting arrays have 1105920 elements and show a peak at 1kHz and another at -1kHz. After ffting a real function, people often double the positive frequency peak and ignore the negative frequency peak, which is all right as long as you keep in mind what the fft is actually doing.

hmax = 1105920
fs = 44100
% 1k puretone in time domain
  delt = 1/fs;
  f0 = 1000;
  T = 5;             % time record (sec)               
  N = T*44100;
  t = (0:N-1)*delt;
  puretone = (1/sqrt(2))*sin(2*pi*f0*t);  
% back to fft code
  fft_puretone = abs(fftshift(fft(puretone)/N));
  Npad = hmax-N;
  fft_puretone = [zeros(1,Npad/2) fft_puretone zeros(1,Npad/2)];
  delf = 1/(N*delt);   % fundamental rule for fft
  freq = (-hmax/2:hmax/2-1)*delf;
  figure(1)
  plot(freq,fft_puretone)

Related Solutions

MATLAB: Matlab crashes when I paste a long segment of code into an m-file

I can confirm this for R2009a/64/Win7. I've pasted it as multiple lines.

I have opened a new M-file in the editor. Then I've opened the same file using WordPad and inserted the code there. After saving Matlab's editor updated its display and Matlab crashed again. Re-opening Matlab did not crash, but leaves the editor window open. The last semi-colon of the pasted code is marked by M-Lint with the warning:

 An M-Lint problem analyzing this file caused your previous MATLAB session to terminate unexpectedly.

(Very interesting. I'm curious where Matlab stores this information!)

Now I've created an M-Lint report for the file and I got 1440 messages about the usage of "i" and "j" as variables. Finally I've replaced them by "ii" and "jj", copied it and pasted the block to a new file - and Matlab crashed again.

Besides the obvious method to upgrade your Matlab, I suggest to use a less brute way of programming. I boldly claim without a even trying it: This function will be very slow.

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

Best Answer

Related Solutions

MATLAB: Matlab crashes when I paste a long segment of code into an m-file

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

Related Question