[help] FFT and spectrogram

[febty febriani]

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/address@hidden/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