mirror of
https://github.com/cookiengineer/audacity
synced 2025-05-08 15:52:53 +02:00
40 lines
1.3 KiB
Matlab
40 lines
1.3 KiB
Matlab
% program to test 2d real fft
|
|
|
|
% let user select file then open it
|
|
[fname, pname] = uigetfile('*.dr2', 'select conv file');
|
|
cd(pname);
|
|
fidout=fopen(fname,'r');
|
|
|
|
% read header info
|
|
N=fread(fidout,1,'long');
|
|
M=fread(fidout,1,'long');
|
|
% read in data
|
|
%status=fseek(fidout,Nheader,'bof');
|
|
a=fread(fidout,N*M,'float');
|
|
a=reshape(a,N,M);
|
|
c=fread(fidout,N*M,'float');
|
|
c=reshape(c,N,M);
|
|
c=c(1:2:N,:)+j*c(2:2:N,:);
|
|
fclose(fidout);
|
|
|
|
c2=fft2(a);
|
|
|
|
% Remember Matlab is column major order, C is row major order
|
|
% Matlab starts with index of 1, f and k start at 0
|
|
|
|
%%%% check the real transforms of DC and Nyquest frequencies
|
|
% check the four real values
|
|
maxerr=abs(real(c(1,1))-c2(1,1)); % 0 f, 0 k
|
|
maxerr=max(maxerr,abs(imag(c(1,1))-c2(1,M/2+1))); % 0 f, M/2 k
|
|
maxerr=max(maxerr,abs(real(c(1,M/2+1))-c2(N/2+1,1)));% N/2 f, 0 k
|
|
maxerr=max(maxerr,abs(imag(c(1,M/2+1))-c2(N/2+1,M/2+1)));% N/2 f, M/2 k
|
|
%check the rest of the pos wavenumbers of DC and Nyquest frequencies
|
|
maxerr=max(maxerr,max(abs(c(1,2:M/2)-c2(1,2:M/2)))); %f = 0, k=1 to M/2-1
|
|
maxerr=max(maxerr,max(abs(c(1,M/2+2:M)-c2(N/2+1,2:M/2))));%f = N/2, k=1 to M/2-1
|
|
%%%%
|
|
|
|
% check all the other positive frequencies at all wavenumbers
|
|
% f from 1 to N/2-1, k from 0 to M-1 (wraps around through negative k)
|
|
maxerr=max(maxerr,max(max(abs(c(2:N/2,:)-c2(2:N/2,:)))))
|
|
|