Why is frequency of graph differing?

Question

Niklas Kurz 2021-1-17

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/719430-why-is-frequency-of-graph-differing

回答： Mathieu NOE 2021-1-19

采纳的回答： Mathieu NOE

在 MATLAB Online 中打开

I like to graph the frequency of 440 Hz and analyse it

I use the following code:

figure(1)
sf = 1000; %sampling
f = 440;
a = 0.5;
t = 0:1/sf:1; 
x = a*sin(2*pi*f*t);
sound(x,sf)
figure(1)
plot(t,x) %amplitude-time graph
axis([0 1/f min(x) max(x)])
figure(2) %Fourier analysis 
n = length(t);
xhat = fft(x,n);
PSD = xhat.*conj(xhat)/n;
freq = 1/(dt*n)*(0:n);
L = 1:floor(n/2);
colormap jet(20)
plot(freq(L),PSD(L),'LineWidth',2.5)

First of: why is it, that the Frequency of amp-time Graph is jagged? If I increase the sampling frequency, why is, that the actual frequency changing?

6 个评论
显示 4更早的评论隐藏 4更早的评论

Mathieu NOE 2021-1-18

在 MATLAB Online 中打开

so I added dt = 1/sf

but the output of the code is fine and the graph is Ok

what's the problem ?

just in case it may interest you , same kina fft analysis + findpeaks demo for notch filter demo

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% dummy data
Fs = 1000;
samples = 5000;
t = (0:samples-1)'*1/Fs;
signal = cos(2*pi*50*t)+cos(2*pi*440*t)+1e-3*rand(samples,1); % two sine + some noise
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = Fs;    % so delta f = 1 Hz
Overlap = 0.75;
w = hanning(NFFT);     % Hanning window / Use the HANN function to get a Hanning window which has the first and last zero-weighted samples.
%% notch filter section %%%%%%
% H(s) = (s^2 + 1) / (s^2 + s/Q + 1)
fc = 50;       % notch freq
wc = 2*pi*fc;
Q = 5;         % adjust Q factor for wider (low Q) / narrower (high Q) notch
                % at f = fc the filter has gain = 0
w0 = 2*pi*fc/Fs;
alpha = sin(w0)/(2*Q);
            b0 =   1;
            b1 =  -2*cos(w0);
            b2 =   1;
            a0 =   1 + alpha;
            a1 =  -2*cos(w0);
            a2 =   1 - alpha;
            
% analog notch (for info)
num1=[1/wc^2 0 1];
den1=[1/wc^2 1/(wc*Q) 1];
% digital notch (for info)
num1z=[b0 b1 b2];
den1z=[a0 a1 a2];
freq = linspace(fc-1,fc+1,200);
[g1,p1] = bode(num1,den1,2*pi*freq);
g1db = 20*log10(g1);
[gd1,pd1] = dbode(num1z,den1z,1/Fs,2*pi*freq);
gd1db = 20*log10(gd1);
figure(1);
plot(freq,g1db,'b',freq,gd1db,'+r');
title(' Notch: H(s) = (s^2 + 1) / (s^2 + s/Q + 1)');
legend('analog','digital');
xlabel('Fréquence (Hz)');
ylabel(' dB')
 
% now let's filter the signal 
signal_filtered = filtfilt(num1z,den1z,signal);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display : averaged FFT spectrum before / after notch filter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[freq,fft_spectrum] = myfft_peak(signal, Fs, NFFT, Overlap);
sensor_spectrum_dB = 20*log10(fft_spectrum);% convert to dB scale (ref = 1)
% demo findpeaks
df = mean(diff(freq));
[pks,locs]= findpeaks(sensor_spectrum_dB,'SortStr','descend','NPeaks',2);
[freq,fft_spectrum_filtered] = myfft_peak(signal_filtered, Fs, NFFT, Overlap);
sensor_spectrum_filtered_dB = 20*log10(fft_spectrum_filtered);% convert to dB scale (ref = 1)
figure(2),semilogx(freq,sensor_spectrum_dB,'b',freq,sensor_spectrum_filtered_dB,'r');grid
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
legend('before notch filter','after notch filter');
xlabel('Frequency (Hz)');ylabel(' dB')
text(freq(locs)+1,pks,num2str(freq(locs)))
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

Mathieu NOE 2021-1-18

在 MATLAB Online 中打开

no need for irritations anywhere

the PSD concept is valid only for random type signals , not sinus

mathematically speaking, a sinus has a PSD of infinite amplitude and zero bandwith, simply because a sine wave has its energy concentrated only at one discrete frequency and not spread accross a wider freq range

if you test your code with random signals, you should see that the peak amplitude stays the same whatever the sampling frequency. But this does not hold if you test now with a sine wave

for a sine wave , the psd concept must be replaced by power or rms or peak amplitude

in summary , a FFT analyser will operate this way :

it will first compute the ftt.*conj(fft) power spectrum

if you ask for "PSD" (knowing that you are looking at random type signals) , it will do : PSD = power / nfft (as you did)

if you ask for "RMS" (knowing that you are looking at sine / multi sine type signals) , it will do : RMS = sqrt(power) and there is no division by nfft

so I modified a bit your code to test the different scenarios :

sf = 2200; %sampling
f = 440;
a = 0.5;
t = 0:1/sf:1; 
% x = a*sin(2*pi*f*t); % to test the PSD = 2*abs(xhat)/n
x = a*randn(size(t)); % to test the PSD = xhat.*conj(xhat)/n
sound(x,sf)
figure(1)
plot(t,x) %amplitude-time graph
axis([0 1/f min(x) max(x)])
figure(2) %Fourier analysis 
n = sf; % so df = sf/n = 1 Hz
xhat = fft(x,n);
PSD = xhat.*conj(xhat)/n; % ok for random type signals
% PSD = 2*abs(xhat)/n; % ok for sine  type signals
freq = 1/(1/sf*n)*(0:n);
L = 1:floor(n/2);
colormap jet(20)
plot(freq(L),PSD(L),'LineWidth',2.5)

Niklas Kurz 2021-1-18

can you post this as answere pls so I can accept .

Mathieu NOE 2021-1-19

will do

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Mathieu NOE 2021-1-19

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/719430-why-is-frequency-of-graph-differing#answer_601034

在 MATLAB Online 中打开

hi - official answer this time

the PSD concept is valid only for random type signals , not sinus

mathematically speaking, a sinus has a PSD of infinite amplitude and zero bandwith, simply because a sine wave has its energy concentrated only at one discrete frequency and not spread accross a wider freq range

if you test your code with random signals, you should see that the peak amplitude stays the same whatever the sampling frequency. But this does not hold if you test now with a sine wave

for a sine wave , the psd concept must be replaced by power or rms or peak amplitude

in summary , a FFT analyser will operate this way :

it will first compute the ftt.*conj(fft) power spectrum

if you ask for "PSD" (knowing that you are looking at random type signals) , it will do : PSD = power / nfft (as you did)

if you ask for "RMS" (knowing that you are looking at sine / multi sine type signals) , it will do : RMS = sqrt(power) and there is no division by nfft

so I modified a bit your code to test the different scenarios :

sf = 2200; %sampling
f = 440;
a = 0.5;
t = 0:1/sf:1; 
% x = a*sin(2*pi*f*t); % to test the PSD = 2*abs(xhat)/n
x = a*randn(size(t)); % to test the PSD = xhat.*conj(xhat)/n
sound(x,sf)
figure(1)
plot(t,x) %amplitude-time graph
axis([0 1/f min(x) max(x)])
figure(2) %Fourier analysis 
n = sf; % so df = sf/n = 1 Hz
xhat = fft(x,n);
PSD = xhat.*conj(xhat)/n; % ok for random type signals
% PSD = 2*abs(xhat)/n; % ok for sine  type signals
freq = 1/(1/sf*n)*(0:n);
L = 1:floor(n/2);
colormap jet(20)
plot(freq(L),PSD(L),'LineWidth',2.5)

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Why is frequency of graph differing?

6 个评论
显示 4更早的评论隐藏 4更早的评论

采纳的回答

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

Community Treasure Hunt

Why is frequency of graph differing?

6 个评论 显示 4更早的评论隐藏 4更早的评论

采纳的回答

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

Community Treasure Hunt

6 个评论
显示 4更早的评论隐藏 4更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论