!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Sergio Steffe' - Laboratorio Sperimentale di Matematica Computazionale
! Dipartimento di Matematica - Universita' di Pisa - AA 2002/2003
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!!!!!!!!!!!!!!! descrizione !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! esempio di filtraggio usando la FFT 
! programma principale e varie subroutine
! i files sonori sono  stereo 16 bit PCM big endian 44100Hz
! i dati del filtro sono nel file dati.filtro
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
subroutine leggiblocco(unita,blocco,posiniz,blcands,blcansn,letti)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!legge un blocco di dati
!
!unita = numero della unita di lettura
!blocco= numero di records da leggere
!posiniz=posizione del primo record
!blcands,blcansn= vettori double precision dei dati del canale ds e sn lunghi blocco
!letti=numero dei record davvero letti
!
!in caso di fine file od errore, vengono letti 0.0D0 nei dati successivi del blocco
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
IMPLICIT NONE
INTEGER 				:: unita,blocco,posiniz,letti
DOUBLE PRECISION,DIMENSION(blocco)	:: blcands,blcansn
INTEGER					:: i,j,ioerr
INTEGER(2)				:: ii,jj,iids,iisn
!
!legge un blocco di dati
!
do i=posiniz,posiniz+blocco
read(unit=unita,rec=i,IOSTAT=ioerr)ii,jj
IF (ioerr /= 0) EXIT
call MVBITS(ii,0,8,iids,8)
call MVBITS(ii,8,8,iids,0)
call MVBITS(jj,0,8,iisn,8)
call MVBITS(jj,8,8,iisn,0)
j=i-posiniz+1
blcands(j)=iids
blcansn(j)=iisn
end do
!se non sos stati letti tutti i dati, metti i rimanenti a 0
do i=j,blocco
blcands(j)=0.0
blcansn(j)=0.0
end do
letti=j
!
end subroutine leggiblocco
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
subroutine scriviblocco(unita,blocco,posiniz,blcands,blcansn,maxrecord)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!unita = numero della unita di scrittura
!blocco= numero di records da scrivere
!posiniz=posizione del primo record
!blcands,blcansn= vettori double precision dei dati del canale ds e sn lunghi blocco
!maxrecord=numero assoluto del massimo record scrivibile
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!
IMPLICIT NONE
INTEGER                                 :: unita,blocco,posiniz,maxrecord
DOUBLE PRECISION,DIMENSION(blocco)      :: blcands,blcansn
INTEGER                                 :: i,j,nmax
INTEGER(2)                              :: ii,iids,iisn
!
!scrive un blocco di dati
!
nmax=min(posiniz+blocco,maxrecord)
!
do i=posiniz,nmax
j=i-posiniz+1
ii=NINT(blcands(j))
call MVBITS(ii,0,8,iids,8)
call MVBITS(ii,8,8,iids,0)
ii=NINT(blcansn(j))
call MVBITS(ii,0,8,iisn,8)
call MVBITS(ii,8,8,iisn,0)
!
write(unit=unita,rec=i) iids,iisn
end do
end subroutine scriviblocco
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
subroutine lunghezza(unita,numrecord)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!unita = numero della unita di scrittura
! si suppone l'unita aperta OLD,UNFORMATTED,DIRECT,RECL=4
!
!numrecord=numero assoluto dei records nel file
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
IMPLICIT NONE
INTEGER                                 :: unita,numrecord
INTEGER 				:: i,ioerr
INTEGER(2)				:: ii,jj
!
numrecord=0
ioerr=0
i=0
!
DO
i=i+1
read(unit=unita,rec=i,IOSTAT=ioerr)ii,jj
IF (ioerr /= 0) EXIT
END DO
!
numrecord=i-1
!
end subroutine lunghezza
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
subroutine settafiltro(nfft,fsample,fftfiltro)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!Il file dati.filtro contiene NDATI coppie di dati
! frequenza  attenuazione
! frequenza da 0.0 a 22050.0
! attenuazione da 0.0 a 1.0
!
!fftfiltro contiene nfft/2 elementi
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
IMPLICIT NONE
INTEGER                                 	:: nfft
DOUBLE PRECISION				:: fsample
DOUBLE PRECISION,DIMENSION(nfft)		:: fftfiltro
INTEGER                                		:: ndati,i,kk,j
DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE	:: fre,smorz 
DOUBLE PRECISION				:: x
!
open(unit=12,file='dati.filtro',status='OLD')
read(12,*) ndati
ALLOCATE(fre(ndati))
ALLOCATE(smorz(ndati))
do i=1,ndati
read(12,*)fre(i),smorz(i)
end do
!
!
!fsample frequenza campionamento
!blocco FFT di nfft dati
!i termini della FFT vanno da -nfft/2 a +nfft/2
!al termine n corrisponde frequenza=(n*fsample)/nfft
!esempio n=nfft/2 ---> fsample/2 = Nyquist
!
do i=1,nfft/2
!
!frequenza corrispondente
x=i
x=fsample*x/nfft
!
kk=1
do j=1,ndati-1
if ((fre(j) < x).and.(x < fre(j+1))) THEN 
 kk=j
end if
end do
!
fftfiltro(i)=smorz(kk)+(smorz(kk+1)-smorz(kk))*(x-fre(kk))/(fre(kk+1)-fre(kk))
end do
!
end subroutine settafiltro
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine copiaxfft(nfft,blcands,blcansn,fftds,fftsn)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 
! copia i blocchi nei buffer per la fft
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
IMPLICIT NONE
INTEGER					:: nfft
DOUBLE PRECISION,DIMENSION(nfft)	:: fftds,fftsn
DOUBLE PRECISION,DIMENSION(nfft/2)   	:: blcands,blcansn
INTEGER					:: i,j
!
do i=1,NFFT/4-1
fftds(i)=0.0D0
fftsn(i)=0.0D0
end do
!
do i=NFFT/4,3*NFFT/8-1
j=1+i-NFFT/4
fftds(i)=(j-1)*blcands(j)/(NFFT/8)
fftsn(i)=(j-1)*blcansn(j)/(NFFT/8)
end do
!
do i=3*NFFT/8,5*NFFT/8-1
j=1+i-NFFT/4
fftds(i)=blcands(j)
fftsn(i)=blcansn(j)
end do
!
do i=5*NFFT/8,3*NFFT/4-1
j=1+i-NFFT/4
fftds(i)=(NFFT/2-j-1)*blcands(j)/(NFFT/8)
fftsn(i)=(NFFT/2-j-1)*blcansn(j)/(NFFT/8)
end do
!
do i=3*NFFT/4,NFFT
fftds(i)=0.0D0
fftsn(i)=0.0D0
end do
!
end subroutine copiaxfft
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!
program filtro
!
!!!!!!!!!!!!!!! descrizione !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!prende un file stereo 16 bit PCM big endian e 
!          lo filtra in base ai parametri del file filtro.dati
!	   che contiene per i vari range di frequenze un fattore
!          di smorzamento:
!	   frequenza-inizio-range frequenza-fine-range fattore (tra 0 ed 1)
!	   i range non citati hanno fattore di default 1
!
!
!!!!!!!!!!!!!!! dichiarazioni !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
IMPLICIT NONE
!
INTEGER,PARAMETER ::  iduemax=32767
!mi tengo sotto il minimo valore 32767 per non avere guai con l'arrotondamento
!
!frequenza di campionamento
DOUBLE PRECISION :: fsample=44100.0
!
!numero campionature
INTEGER :: ntot
!
!posizione iniziale lettura e scrittura
INTEGER :: posin, posout
!
!interi usati per contatori
INTEGER i,j
!
!lunghezza blocchi da elaborare
INTEGER :: BLOCCO
!
!vettori dei dati letti a blocco
DOUBLE PRECISION, DIMENSION(:),ALLOCATABLE :: blcands,blcansn,pezzods,pezzosn
!
!tempo totale
DOUBLE PRECISION  ttot
!
!dati letti
INTEGER letti
!
!nomi del file di  output
CHARACTER(len=25) FILENAMEOUT
!
!posfisso da aggiungere al nome file di output
CHARACTER(len=4) :: POSTF='.cdr'
!
!variabile per lettura nome file 
CHARACTER(len=20) FILENAME1
!
!nome del file di input
CHARACTER(len=64) FILENAMEIN
!
!segnale sonoro a ds e sn 
!DOUBLE PRECISION soundds,soundsn
!
!per la trasformata di fourier
DOUBLE PRECISION, DIMENSION(:),ALLOCATABLE :: fftds,fftsn,fftmp,fftfiltro
INTEGER NFFT
!
!

!!!!!!!!!!!!!! inizio !!!!!!!!!!!!!!!!!!!!!!!!!!!
!si aprono i files
!REPLACE crea o cancella se esistente
!UNFORMATTED usa formato interno per salvare i dati
!DIRECT usa l'accesso ad un qualsiasi record
!RECL e' la lunghezza di ogni record (4 bytes)
!
!
write(*,*)'Nome del file di input:'
read(*,*)FILENAMEIN
!write(*,*)'hai detto: ',FILENAMEIN
!
open(unit=11,file=FILENAMEIN,status='OLD',form='UNFORMATTED',access='DIRECT',RECL=4,ERR=101)
!
write(*,*)'Nome del file di output (max 20 char senza spazi): '
read(*,*) FILENAME1
FILENAMEOUT = FILENAME1(1:SCAN( FILENAME1,' ',.FALSE.)-1)//POSTF
!write(*,*)'hai detto: ',FILENAMEOUT
!
open(unit=10,file=FILENAMEOUT,status='REPLACE',form='UNFORMATTED',access='DIRECT',RECL=4)
!
!!!!!!!!!!!!!!!!!!!!allocazioni matrici!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
BLOCCO=2**14
ALLOCATE(blcands(BLOCCO))
ALLOCATE(blcansn(BLOCCO))
NFFT=BLOCCO*2
ALLOCATE(fftds(NFFT))
ALLOCATE(fftsn(NFFT))
ALLOCATE(fftmp(2*NFFT+15))
ALLOCATE(fftfiltro(NFFT/2))
ALLOCATE(pezzods(BLOCCO/4))
ALLOCATE(pezzosn(BLOCCO/4))
!
!
call settafiltro(NFFT,fsample,fftfiltro)
!
!!!!!!!!!!!!!!!!!!inizio operazioni calcolo!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!calcola ntot
write(*,*)'Conto il numero dei records di ',FILENAMEIN
call lunghezza(11,ntot)
write (*,*) 'record letti =',ntot
!
!puntatore assoluto lettura e scrittura
!
!inizializza FFT
call dffti(NFFT,fftmp)
!
posin=1
posout=1
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! leggo un blocco di dati
!
letti=0
call leggiblocco(11,blocco,posin,blcands,blcansn,letti)
posin=posin+letti
!
!salvo pezzi iniziali scalati in modo opposto alla fft
do i=1,blocco/4
pezzods(i)=(1.0 - i*4.0/blocco)*blcands(i)
pezzosn(i)=(1.0 - i*4.0/blocco)*blcansn(i)
end do
!
call copiaxfft(nfft,blcands,blcansn,fftds,fftsn)
!
call dfftf(NFFT,fftds,fftmp)
call dfftf(NFFT,fftsn,fftmp)
!
!filtraggio qui
fftds(1)=fftds(1)*fftfiltro(1)
fftsn(1)=fftsn(1)*fftfiltro(1)
fftds(nfft)=fftds(nfft)*fftfiltro(nfft/2)
fftsn(nfft)=fftsn(nfft)*fftfiltro(nfft/2)
do i=1,nfft/2-1
fftds(2*i)=fftds(2*i)*fftfiltro(i)
fftds(2*i+1)=fftds(2*i+1)*fftfiltro(i)
fftsn(2*i)=fftsn(2*i)*fftfiltro(i)
fftsn(2*i+1)=fftsn(2*i+1)*fftfiltro(i)
end do
!
call dfftb(NFFT,fftds,fftmp)
call dfftb(NFFT,fftsn,fftmp)
!
!recupera blocco dopo fft e rinormalizza
!
do i=1,blocco
j=i+NFFT/4
blcands(i)=fftds(j)/(NFFT)
blcansn(i)=fftsn(j)/(NFFT)
end do
!
! il primo pezzetto del primo blocco va unito a quello non sfumato non filtrato
do i=1,blocco/4
blcands(i)=blcands(i)+pezzods(i)
blcansn(i)=blcansn(i)+pezzosn(i)
end do
!
!va salvato ora il secondo pezzetto del primo blocco
do i=1,blocco/4
pezzods(i)=blcands(i+(3*blocco)/4)
pezzosn(i)=blcansn(i+(3*blocco)/4)
end do
!
!si scrive ora la parte fatta
!
call scriviblocco(10,3*blocco/4,posout,blcands,blcansn,ntot)
posout=posout+3*blocco/4
!
!fine del primo passo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!passi da iterare fino a fine file si legge due volte un pezzetto, pazienza
DO while (posin +3*blocco/4 < ntot)
letti=0
call leggiblocco(11,blocco,posin-blocco/4,blcands,blcansn,letti)
posin=posin-blocco/4+letti
call copiaxfft(nfft,blcands,blcansn,fftds,fftsn)
call dfftf(NFFT,fftds,fftmp)
call dfftf(NFFT,fftsn,fftmp)
!filtraggio qui
fftds(1)=fftds(1)*fftfiltro(1)
fftsn(1)=fftsn(1)*fftfiltro(1)
fftds(nfft)=fftds(nfft)*fftfiltro(nfft/2)
fftsn(nfft)=fftsn(nfft)*fftfiltro(nfft/2)
do i=1,nfft/2-1
fftds(2*i)=fftds(2*i)*fftfiltro(i)
fftds(2*i+1)=fftds(2*i+1)*fftfiltro(i)
fftsn(2*i)=fftsn(2*i)*fftfiltro(i)
fftsn(2*i+1)=fftsn(2*i+1)*fftfiltro(i)
end do
!
call dfftb(NFFT,fftds,fftmp)
call dfftb(NFFT,fftsn,fftmp)
!recupera blocco dopo fft e rinormalizza
do i=1,blocco
j=i+NFFT/4
blcands(i)=fftds(j)/(NFFT)
blcansn(i)=fftsn(j)/(NFFT)
end do
! il primo pezzetto  blocco va unito a quello finale precedente
do i=1,blocco/4
blcands(i)=blcands(i)+pezzods(i)
blcansn(i)=blcansn(i)+pezzosn(i)
end do
!
!va salvato ora il pezzetto finale del blocco
do i=1,blocco/4
pezzods(i)=blcands(i+(3*blocco)/4)
pezzosn(i)=blcansn(i+(3*blocco)/4)
end do
!
!si scrive ora la parte fatta
!
call scriviblocco(10,3*blocco/4,posout,blcands,blcansn,ntot)
posout=posout+3*blocco/4
!
!fine del while
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!ultimo passaggio ora posin +3*blocco/4 > ntot
!
letti=0
call leggiblocco(11,blocco,posin-blocco/4,blcands,blcansn,letti)
posin=posin-blocco/4+letti
call copiaxfft(nfft,blcands,blcansn,fftds,fftsn)
call dfftf(NFFT,fftds,fftmp)
call dfftf(NFFT,fftsn,fftmp)
!filtraggio qui
fftds(1)=fftds(1)*fftfiltro(1)
fftsn(1)=fftsn(1)*fftfiltro(1)
fftds(nfft)=fftds(nfft)*fftfiltro(nfft/2)
fftsn(nfft)=fftsn(nfft)*fftfiltro(nfft/2)
do i=1,nfft/2-1
fftds(2*i)=fftds(2*i)*fftfiltro(i)
fftds(2*i+1)=fftds(2*i+1)*fftfiltro(i)
fftsn(2*i)=fftsn(2*i)*fftfiltro(i)
fftsn(2*i+1)=fftsn(2*i+1)*fftfiltro(i)
end do
!
call dfftb(NFFT,fftds,fftmp)
call dfftb(NFFT,fftsn,fftmp)
!recupera blocco dopo fft e rinormalizza
do i=1,blocco
j=i+NFFT/4
blcands(i)=fftds(j)/(NFFT)
blcansn(i)=fftsn(j)/(NFFT)
end do
! il primo pezzetto  blocco va unito a quello finale precedente
do i=1,blocco/4
blcands(i)=blcands(i)+pezzods(i)
blcansn(i)=blcansn(i)+pezzosn(i)
end do
!
!
!si scrive ora la parte finale
!
call scriviblocco(10,blocco,posout,blcands,blcansn,ntot)
!
!
close(unit=10)
close(unit=11)
!
!
ttot=ntot
ttot=ttot/fsample
write(*,501)ttot,FILENAMEOUT
501 FORMAT('Scritti',F5.1,' secondi suono (PCM 2x16 bit big endian) nel file: ',A)
GOTO 202
!
101 write(*,*)'ERRORE aprendo in lettura il file: ',FILENAMEIN
GOTO 202
!
102  write(*,*)'ERRORE leggendo il file: ',FILENAMEIN
ttot=i
ttot=ttot/fsample
write(*,*)'(dopo ',ttot,' secondi di audio)'
GOTO 202
!
105  write(*,*)'ERRORE: aprendo il file dei dati del filtro '
GOTO 202
!
202 CONTINUE
end program filtro
