Get spectrum from autocorrelation function using fft()

22 visualizaciones (últimos 30 días)
Ivan Volodin
Ivan Volodin el 22 de Mzo. de 2017
Editada: Ivan Volodin el 22 de Mzo. de 2017
Hello!
I would like to observe Fourier Spectrum, got through Autocorrelation function. This is my code:
clear;clc;
%%generate x-axe
first_step=0;
step_t=0.01;
last_step=1;
t=first_step:step_t:last_step;
%%signal and fourier transform
y= 2*sin(2*pi*60*t) ;
figure(1)
plot(t,y);title('signal')
fourier=fft(y );
N_=length(fourier);
f_ = (1 / step_t) * ( 0: (N_/2) ) / N_;
a=(fourier.*conj(fourier))/(N_*N_); % Spectral Density get half
a = a(1:N_ /2+1);
a(2:end-1) = 4*a(2:end-1);
figure(2);
plot(f_,a);title('Spectrum Power through signal');
%%through AutoCorrelation function
tau=0.01;
y1= 2*sin(2*pi*60*(t+tau)) ;
Y=y1.*y;
AutoCorr_func = cumtrapz(Y,t);
figure(3);
plot(t,AutoCorr_func); title('AutoCorrelation function');
fourier=fft(AutoCorr_func);
N_=length(fourier);
f_ = (1 / step_t) * ( 0: (N_/2) ) / N_;
a=fourier.*conj(fourier) ; % Spectral Density
a = a(1:N_ /2+1); % spectrum is even, so get half
a(2:end-1) = a(2:end-1); % multiply by 4, since spectrum is square of amplitude, and we have two even halves
figure(4)
plot(f_ ,a );
title('Spectrum through AutoCorrelation funciton');
and I get this
so a few questions:
1) why I get so strang ACF? It has to be periodic, since y and y1 are strongly periodic. Is sonething wrong with it? maybe here
Y=y1.*y;
AutoCorr_func = cumtrapz(Y,t);
but Acf is an integral and I am trying to get an array numerically.
2) My tries to get power density through Acf are the same as in the simple case, without using Acf, but Figure 2 gives correct resultL frequency = 42 and power is square of amplitude, in my case it should be 4, but I have around 2. What is the problem...?
3) How to get correct power density through Acf, what am I doing wrong?
Any help would be greatly appreciate! Thank you in advance!

Respuestas (0)

Categorías

Más información sobre Spectral Measurements en Help Center y File Exchange.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by