Prev: problem assigning variables with eval - nameclash with existing functions
Next: vim as editor for .m files
From: Agus Nugroho on 1 Jul 2010 10:01
Dear comrade in arms,

I would like to ask you concerning the problem that I am facing now and I hope you do not mind to help me.

I have tested a sinusoidal signal to find it's auto spectrum with algorithms (X- axis represents frequency and Y-axis represents Amplitude). The X-axis (frequency) of the curve is correct, however, I am not really sure about the Y-Axis (amplitude). Could you help me by seeing my algorithm below whether that is correct to find its amplitude? :


-------------------------------------------------------------------------
Fs=20;
Ts=1/Fs;
t=0:Ts:12.75;
n=length(t);
y=sin(2*pi*0.3515*t); %generate sinusoidal wave with freq 0.3515
figure(1)
plot(t,y)

[YfreqD,freqRng]=positiveFFT(y,Fs);

figure(2)
YfreqDb=2*abs(YfreqD); % amplitude scale

Auto_sp=2*(conj(YfreqDb).*YfreqDb); %auto spectrum signal

Auto_scale=sqrt(2*Auto_sp)*(2/256); %Auto spectrum scale, 256 is nfft

plot(freqRng, Auto_scale)

-----------------------------------------------------------------------------------

the algorithm for function :

function[YfreqD,freqRng]=positiveFFT(windowedSignal,Fs)
N=length(windowedSignal);
k=0:N-1;
T=N/Fs;
freqRng=k/T;
YfreqD=fft(windowedSignal)/N;
cutOff=ceil(N/2);
YfreqD=YfreqD(1:cutOff);
freqRng=freqRng(1:cutOff);
-------------------------------------------------------------------

Thank you very much for your help.

Kindest regards,

Agus
From: Wayne King on 1 Jul 2010 11:37
"Agus Nugroho" <agus_nugroho01(a)yahoo.com> wrote in message <i0i73f$m8i$1(a)fred.mathworks.com>...
> Dear comrade in arms,
>
> I would like to ask you concerning the problem that I am facing now and I hope you do not mind to help me.
>
> I have tested a sinusoidal signal to find it's auto spectrum with algorithms (X- axis represents frequency and Y-axis represents Amplitude). The X-axis (frequency) of the curve is correct, however, I am not really sure about the Y-Axis (amplitude). Could you help me by seeing my algorithm below whether that is correct to find its amplitude? :
>
>
> -------------------------------------------------------------------------
> Fs=20;
> Ts=1/Fs;
> t=0:Ts:12.75;
> n=length(t);
> y=sin(2*pi*0.3515*t); %generate sinusoidal wave with freq 0.3515
> figure(1)
> plot(t,y)
>
> [YfreqD,freqRng]=positiveFFT(y,Fs);
>
> figure(2)
> YfreqDb=2*abs(YfreqD); % amplitude scale
>
> Auto_sp=2*(conj(YfreqDb).*YfreqDb); %auto spectrum signal
>
> Auto_scale=sqrt(2*Auto_sp)*(2/256); %Auto spectrum scale, 256 is nfft
>
> plot(freqRng, Auto_scale)
>
> -----------------------------------------------------------------------------------
>
> the algorithm for function :
>
> function[YfreqD,freqRng]=positiveFFT(windowedSignal,Fs)
> N=length(windowedSignal);
> k=0:N-1;
> T=N/Fs;
> freqRng=k/T;
> YfreqD=fft(windowedSignal)/N;
> cutOff=ceil(N/2);
> YfreqD=YfreqD(1:cutOff);
> freqRng=freqRng(1:cutOff);
> -------------------------------------------------------------------
>
> Thank you very much for your help.
>
> Kindest regards,
>
> Agus

Hi Agus, Do you have the Signal Processing Toolbox? If so how about using a spectral analysis object and msspectrum for this?

Ts = 1/20;
t = 0:Ts:12.75;
y = sin(2*pi*5*t);
H = spectrum.periodogram('Hamming');
Hmss = msspectrum(H,y,'Fs',20,'SpectrumType','twosided','CenterDC',true);
plot(Hmss);

Note the power in each of the complex exponentials is
20*log10(1/2) = -6.0206 dB

Just what you would expect it to be since each of the complex exponentials has amplitude 1/2.

Or is there some reason you are trying to do this from scratch? If so write back and I'll show you what you've done wrong. Just don't want you reinventing the wheel so to speak.

Wayne
From: Agus Nugroho on 1 Jul 2010 12:20
Dear Dr. Wayne,

Thank you very much for your valuable help. Concerning my reason to find the auto spectrum from the strech is due to the advice of my supervisor. Our research is to find the evolution of the auto spectrum from time to time, and my supervisor though the best way to investigate this is by doing the methods that I have described previously.

So I tried to tested the methodology from strech by using a sample signal (I found the sample signal from one paper, although they just provide the signal sample and the graph of auto spectrum).

From tested the frequency and its curve shape are correct, but the amplitude is different. Therefore I would like to ask if I made some mistake in my algorithm or not. Perhaps with this, I also could to compare with the Hmss methodology. Your help is highly appreciated.

Thank you very much for your help

Kindest regards,

Agus
From: Wayne King on 1 Jul 2010 12:58
"Agus Nugroho" <agus_nugroho01(a)yahoo.com> wrote in message <i0if87$ohm$1(a)fred.mathworks.com>...
> Dear Dr. Wayne,
>
> Thank you very much for your valuable help. Concerning my reason to find the auto spectrum from the strech is due to the advice of my supervisor. Our research is to find the evolution of the auto spectrum from time to time, and my supervisor though the best way to investigate this is by doing the methods that I have described previously.
>
> So I tried to tested the methodology from strech by using a sample signal (I found the sample signal from one paper, although they just provide the signal sample and the graph of auto spectrum).
>
> From tested the frequency and its curve shape are correct, but the amplitude is different. Therefore I would like to ask if I made some mistake in my algorithm or not. Perhaps with this, I also could to compare with the Hmss methodology. Your help is highly appreciated.
>
> Thank you very much for your help
>
> Kindest regards,
>
> Agus

Hi Agus, I'm a little confused on what your use case is. If you are trying to estimate the PSD--the auto spectrum of a WSS process, then just use the code I gave you in the other post. Here's another example:

Fs = 20;
Ts = 1/Fs;
t = 0:Ts:12.75;
N = length(t);
y = sin(2*pi*5*t);
[Pxx,F] = periodogram(y,[],length(y),20);
subplot(211);
plot(F,Pxx)
% compare this to
subplot(212);
yDFT = 2/(length(y)*20)*abs(fft(y)).^2;
plot(F,yDFT(1:129))

Note the agreement. My other post was actually more detailed and technically correct since here I'm ignoring what I told you about DC and the Nyquist. (again see my other post to you)

If you want to use Fourier methods to estimate the amplitude of a sine wave, that is a different application.

Wayne
From: Agus Nugroho on 1 Jul 2010 16:24
Dear Dr Wayne,

After consultate with my supervisor again, he is agree to use the previous method which is Periodogram. And when I shown him about your methods, he is totally agree. Thank you very much for your help Dr.Wayne, your explanation in the las post also was very clear and detail.

Thank you very much

Kindest regards,

Agus
 |  Next  |  Last
Pages: 1 2
Prev: problem assigning variables with eval - nameclash with existing functions
Next: vim as editor for .m files