[help] FFT & Spectrogram

[febty]

Hi everyone,

Maybe my problem is a trivial thing but I have no idea how to solve it. I must display the FFT time series data for many years in segments. My data contain of maximum five channel. In this email, I just attached the script for channel one. I don't know how Octave can handle my result, so I used gnuplot to display my image (the first script: FFT script for channel one).
I tried to attach my figure, but it was too big and my email was bounced so I put my figure in this link : http://www.flickr.com/photos/40025007@N05/3680231469/. It is OK for me. But when I want to display spectrogram's result, I am not sure GnuPlot can do it. So I change my script like the second script using Octave only (FFT&spectrogram script). I want to display the result like my before result using GnuPlot. But, I didn't get anything, even the error message.


Any kind of help is appreciated.
Thanks very much.

##############
1. FFT script

#!/bin/bash

for year in 2009
do
for month in 01
do
for day in 01
do
for num in 00 01 02 03 04
do

#FFT channel 1

awk 'NR>=1800*'$num'+1 && NR<=1800*'$num'+4096{print $2}' ${year}${month}${day}.1Hz.plr > result.dat

octave <<EOF
fs=1;                          
fh=fopen("result.dat","r");
x=fscanf(fh,"%lf");
x0=hanning(4096);                  
a=(x(1:4096)).*(x0);                  
b=fft(a,8192);                      
c=abs(b(1:4096));                  
c=c/max(max(c));
save -ascii ${year}${month}${day}.1Hz.${num}fft1.dat c  
d=((1:4096)/4096)*(fs/2);              
e=d';                          
save -ascii ${year}${month}${day}.1Hz.${num}trans1.dat e  
EOF

paste ${year}${month}${day}.1Hz.${num}trans1.dat ${year}${month}${day}.1Hz.${num}fft1.dat > ${year}${month}${day}.1Hz.${num}compile1.dat

gnuplot <<EOF

set term postscript enhanced color
set output "electric-fft-${year}${month}${day}a.1Hz.ps"

set logscale y
set grid
set ytics out font "Times New Roman, 8"
set ylabel "Amplitude" 2,0 font "Times New Roman, 10"
set format y '10^{%L}'
set xlabel "Frekuensi (Hz)" 0,1 font "Times New Roman, 10"
set xrange [0:0.5]
set xtics 0, 0.1, 0.5 font "Times New Roman, 8"

set multiplot

#channel 1

set origin 0.0,0.79
set size 0.275,0.21
set title "E1-00" 0,-1 font "Times New Roman, 10"
plot "${year}${month}${day}.1Hz.00compile1.dat" using 1:2 with lines notitle "${year}${month}${day}.1Hz.00compile1.dat"

set origin 0.0, 0.59
set size 0.275,0.21
set title "E1-01" 0,-1 font "Times New Roman, 10"
plot "${year}${month}${day}.1Hz.01compile1.dat" using 1:2 with lines notitle "${year}${month}${day}.1Hz.01compile1.dat"

set origin 0.0, 0.39
set size 0.275,0.21
set title "E1-02" 0,-1 font "Times New Roman, 10"
plot "${year}${month}${day}.1Hz.02compile1.dat" using 1:2 with lines notitle "${year}${month}${day}.1Hz.02compile1.dat"

set origin 0.0, 0.19
set size 0.275,0.21
set title "E1-03" 0,-1 font "Times New Roman, 10"
plot "${year}${month}${day}.1Hz.03compile1.dat" using 1:2 with lines notitle "${year}${month}${day}.1Hz.03compile1.dat"

set origin 0.0, 0.0
set size 0.275,0.21
set title "E1-04" 0,-1 font "Times New Roman, 10"
plot "${year}${month}${day}.1Hz.04compile1.dat" using 1:2 with lines notitle "${year}${month}${day}.1Hz.04compile1.dat"

unset multiplot
EOF

done
done
done
done

#################

2. FFT&spectrogram script

#!/bin/bash

for year in 2009
do
for month in 01
do
for day in 01
do
for num in 00 01 02 03 04
do

#FFT channel 1

awk 'NR>=1800*'$num'+1 && NR<=1800*'$num'+4096{print $2}' ${year}${month}${day}.1Hz.plr > result.dat

octave <<EOF
Fs=1;                          
fh=fopen("result.dat","r");
x=fscanf(fh,"%lf");
x0=hanning(4096);                  
a=(x(1:4096)).*(x0);                  
b=fft(a,8192);                      
c=abs(b(1:4096));
c=c/max(max(c))    ;                  
d=((1:4096)/4096)*(Fs/2);                                      
subplot(2,1,1);plot(d,c);
subplot(2,2,1);imagesc(flipud(20*log10(c)));
EOF

done
done
done
done

--

******
febty febriani
Indonesian Institute of Sciences
Research Center for Physics
Kompleks PUSPIPTEK
Serpong, Indonesia, 15314

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Judd Storrs]

On Thu, Jul 2, 2009 at 2:01 AM, febty febriani <febty82@...> wrote:

octave <<EOF
Fs=1;                          
fh=fopen("result.dat","r");
x=fscanf(fh,"%lf");
x0=hanning(4096);                  
a=(x(1:4096)).*(x0);                  
b=fft(a,8192);                      
c=abs(b(1:4096));
c=c/max(max(c))    ;                  
d=((1:4096)/4096)*(Fs/2);                                      
subplot(2,1,1);plot(d,c);
subplot(2,2,1);imagesc(flipud(20*log10(c)));
EOF

What I suspect is happening here is when octave reaches the EOF it exits and doesn't draw the plots. Since there also aren't any commands to save the plots to a file that is why you don't see any output. You can perhaps add a "print" statement to get the plots to file http://www.network-theory.co.uk/docs/octave3/octave_163.html

Or if you would rather watch the plots you can add

drawnow ;
sleep(10) ;

before the EOF to pause 10 seconds before exiting. For the second subplot you may want "subplot(2,1,2)" instead of "subplot(2,2,1)"

Hope this helps,

--judd

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[febty]

2009/7/2 Judd Storrs <storrsjm@...>

On Thu, Jul 2, 2009 at 2:01 AM, febty febriani <febty82@...> wrote:

octave <<EOF
Fs=1;                          
fh=fopen("result.dat","r");
x=fscanf(fh,"%lf");
x0=hanning(4096);                  
a=(x(1:4096)).*(x0);                  
b=fft(a,8192);                      
c=abs(b(1:4096));
c=c/max(max(c))    ;                  
d=((1:4096)/4096)*(Fs/2);                                      
subplot(2,1,1);plot(d,c);
subplot(2,2,1);imagesc(flipud(20*log10(c)));
EOF

What I suspect is happening here is when octave reaches the EOF it exits and doesn't draw the plots. Since there also aren't any commands to save the plots to a file that is why you don't see any output. You can perhaps add a "print" statement to get the plots to file http://www.network-theory.co.uk/docs/octave3/octave_163.html

Or if you would rather watch the plots you can add

drawnow ;
sleep(10) ;

before the EOF to pause 10 seconds before exiting. For the second subplot you may want "subplot(2,1,2)" instead of "subplot(2,2,1)"

Hope this helps,

--judd

thanks very much, it is very helpful for me. I fixed the script. I can see the FFT result. For the spectrogram, I just saw blank image. From Octave prompt, I got error like below :

GNUPLOT (plot_image):  Image grid must be at least 2 x 2.


GNUPLOT (plot_image):  Image grid must be at least 2 x 2.


my script is

###
octave <<EOF
Fs=1;   
n=4096;                       
fh=fopen("result.dat","r");
y=fscanf(fh,"%lf");
window=hanning(n);                   
y1=(y(1:n)).*(window);                   
y2=fft(y1,n*2);   
[S, f, t] = specgram(y, n, Fs, window, length(window)/2);                   
S=abs(y2(1:n));
S=S/max(max(S));                   
x=((1:n)/n)*(Fs/2);                                           
subplot(2,1,1);plot(x,S);axis([0 0.01]);title("200901");xlabel("Frequency");ylabel("Amplitude");
subplot(2,1,2);imagesc(flipud(20*log10(S)));
print -dpng output.png;
drawnow;
sleep(10);
EOF
###

any idea why I got such error?
any help is appreciated. Thanks very much.
--

******
febty febriani
Indonesian Institute of Sciences
Research Center for Physics
Kompleks PUSPIPTEK
Serpong, Indonesia, 15314

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Jaroslav Hajek-2]

On Thu, Jul 2, 2009 at 8:53 AM, Judd Storrs<storrsjm@...> wrote:

> On Thu, Jul 2, 2009 at 2:01 AM, febty febriani <febty82@...> wrote:
>>
>> octave <<EOF
>> Fs=1;
>> fh=fopen("result.dat","r");
>> x=fscanf(fh,"%lf");
>> x0=hanning(4096);
>> a=(x(1:4096)).*(x0);
>> b=fft(a,8192);
>> c=abs(b(1:4096));
>> c=c/max(max(c))    ;
>> d=((1:4096)/4096)*(Fs/2);
>> subplot(2,1,1);plot(d,c);
>> subplot(2,2,1);imagesc(flipud(20*log10(c)));
>> EOF
>
> What I suspect is happening here is when octave reaches the EOF it exits and
> doesn't draw the plots. Since there also aren't any commands to save the
> plots to a file that is why you don't see any output. You can perhaps add a
> "print" statement to get the plots to file
> http://www.network-theory.co.uk/docs/octave3/octave_163.html
>
> Or if you would rather watch the plots you can add
>
> drawnow ;
> sleep(10) ;
>
> before the EOF to pause 10 seconds before exiting. For the second subplot
> you may want "subplot(2,1,2)" instead of "subplot(2,2,1)"
>
> Hope this helps,
>
> --judd

As a side note, I think calling `drawnow' is redundant, because
`sleep' does it intrinsically.

--
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[febty]

2009/7/2 Jaroslav Hajek <highegg@...>

On Thu, Jul 2, 2009 at 8:53 AM, Judd Storrs<storrsjm@...> wrote:
> On Thu, Jul 2, 2009 at 2:01 AM, febty febriani <febty82@...> wrote:
>>
>> octave <<EOF
>> Fs=1;
>> fh=fopen("result.dat","r");
>> x=fscanf(fh,"%lf");
>> x0=hanning(4096);
>> a=(x(1:4096)).*(x0);
>> b=fft(a,8192);
>> c=abs(b(1:4096));
>> c=c/max(max(c))    ;
>> d=((1:4096)/4096)*(Fs/2);
>> subplot(2,1,1);plot(d,c);
>> subplot(2,2,1);imagesc(flipud(20*log10(c)));
>> EOF
>
> What I suspect is happening here is when octave reaches the EOF it exits and
> doesn't draw the plots. Since there also aren't any commands to save the
> plots to a file that is why you don't see any output. You can perhaps add a
> "print" statement to get the plots to file
> http://www.network-theory.co.uk/docs/octave3/octave_163.html
>
> Or if you would rather watch the plots you can add
>
> drawnow ;
> sleep(10) ;
>
> before the EOF to pause 10 seconds before exiting. For the second subplot
> you may want "subplot(2,1,2)" instead of "subplot(2,2,1)"
>
> Hope this helps,
>
> --judd

As a side note, I think calling `drawnow' is redundant, because
`sleep' does it intrinsically.

--
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz

thanks alot. you are right. it was redundant. I eliminated it. but I had no idea for error message "GNUPLOT (plot_image):  Image grid must be at least 2 x 2." for my script above.

any help is appreciated, please.


--

******
febty febriani
Indonesian Institute of Sciences
Research Center for Physics
Kompleks PUSPIPTEK
Serpong, Indonesia, 15314

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Judd Storrs]

On Thu, Jul 2, 2009 at 4:56 AM, febty febriani <febty82@...> wrote:

###
octave <<EOF
Fs=1;   
n=4096;                       

fh=fopen("result.dat","r");

y=fscanf(fh,"%lf");
window=hanning(n);                   
y1=(y(1:n)).*(window);                   
y2=fft(y1,n*2);   
[S, f, t] = specgram(y, n, Fs, window, length(window)/2);                   
S=abs(y2(1:n));
S=S/max(max(S));                   
x=((1:n)/n)*(Fs/2);                                           
subplot(2,1,1);plot(x,S);axis([0 0.01]);title("200901");xlabel("Frequency");ylabel("Amplitude");
subplot(2,1,2);imagesc(flipud(20*log10(S)));
print -dpng output.png;
drawnow;
sleep(10);
EOF
###

I'm not familiar with specgram() but I suspect S is a vector, i.e. 1xN or Nx1 so the imagesc command complains because S isn't an MxN matrix with both N>1 and M>1 (i.e. it's a vector not an image(=matrix)).

You can try adding a disp(size(S)) right before the imagesc call to verify the dimensions.

If you want to display S as an image anyway, I might work to duplicate S to create a 2xN or Nx2 matrix. i.e. something like

S = repmat(S,[1,2]); # if S is Nx1
S = repmat(S,[2,1]); # if S is 1xN

--judd

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave