1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
|
wavefile = 'Testdatei 400 Hz 1.wav';
M = wavread(wavefile); % load data from wave file
x = M(:,1); % wave samples as vector
fs = 16000; % sample frequency
sec = 0.025; % time span
x = x(1:(sec*fs)+1); % samples of time span
wavwrite(x,16000,16,"x.wav")
N = length(x); % number of samples
t = [0:1/fs:sec]; % time vector
nfft = 2^(nextpow2(N)); % Use next highest power of 2 greater than or equal to length(x) to calculate FFT.
NumUniquePts = ceil((nfft+1)/2); % Calculate the numberof unique points
fftx = fft(x,nfft); % Take fft, padding with zeros so that length(fftx) is equal to nfft
save fftfile fftx; % Write fft complex data to file
p = unwrap(angle(fftx)); % phase
ifftx = ifft(fftx);
ifftx = ifftx(1:N);
wavwrite(ifftx,16000,16,"ifftx.wav")
fftx = fftx(1:NumUniquePts); % FFT is symmetric, throw away second half
magX = abs(fftx)/N; % Take the magnitude of fft of x and scale the fft so that it is not a function of the length of x
magX = magX.^2; % Take the square of the magnitude of fft of x.
% Since we dropped half the FFT, we multiply mx by 2 to keep the same energy.
% The DC component and Nyquist component, if it exists, are unique and should not be multiplied by 2.
if rem(nfft, 2) % odd nfft excludes Nyquist point
magX(2:end) = magX(2:end)*2;
else
magX(2:end -1) = magX(2:end -1)*2;
end
f = (0:NumUniquePts-1)*fs/nfft; % This is an evenly spaced frequency vector with NumUniquePts points.
figure(1);
subplot(5,1,1);
plot(t,x);
title('time domain')
xlabel('time');
ylabel('Amplitude');
subplot(5,1,3);
plot(f,magX);
xlabel('f/Hz');
ylabel('Amplitudenspektrum');
title('frequency domain');
subplot(5,1,5);
px = (1:1:length(p));
plot(px,p);
title('phase');
% print -Ppng 'plot.png'
|
