[help] FFT & 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/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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: https://www-old.cae.wisc.edu/pipermail/help-octave/attachments/20090702/ab650bfa/attachment.html