chirpwav file can be dowload from this website httpwwwgrsite

chirp.wav file can be dowload from this website

http://www.grsites.com/archive/sounds/category/10/

for the wav file the frequency rate is 8192.

2. How to use fftshift to shift zero-frequency component to center of spectrum

3. How to calculate the magnitude and phase of DFT components. Plot the double-side amplitude and phase spectrum with respect to frequency (Use Nyquist criteria to define frequency range).

4. How to store the frequency, amplitude and phase parameters in ASCII file (text, excel, csv, etc.)

Hint: audioread, audiowrite, sound, fft, fftshift, stem, save, load, subplot.

Number of Samples

The sampled time waveform input to an FFT determines the computed spectrum. If an arbitrary signal is sampled at a rate equal to f_s over an acquisition time T, N samples are acquired. Compute T with the following equation:

T=N/fs

where

T is the acquisition time

N is the number of samples acquired

f_s is the sampling frequency

Maximum Resolvable Frequency

The sampling rate of a time waveform determines the maximum resolvable frequency. According to the Shannon Sampling Theorem, the maximum resolvable frequency must be half the sampling frequency. To calculate the maximum resolvable frequency, use the following equation:

fmax=fNyquist=fs/2

where

f_max is the maximum resolvable frequency

f_Nyquist is the Nyquist frequency

f_s is the sampling frequency

Solution

clc;
clear all; close all;
%%% sampling frequency %%
Fs=8192;
%%% audio file reading %%%
[x L]= wavread(\'birds001\');
figure
plot(x)

%%% DFT using FFT %%%%
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(x,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);

% Plot single-sided amplitude spectrum.
figure
plot(f,2*abs(Y(1:NFFT/2+1)))
title(\'Single-Sided Amplitude Spectrum of y(t)\')
xlabel(\'Frequency (Hz)\')
ylabel(\'|Y(f)|\')

%%%%% fftshift to shift zero-frequency component to center of spectrum%%%%
shift_centre = fftshift(Y);

%%% magnitude and phase of DFT components %%%
m = abs(Y);
p = angle(Y);
p=rad2deg(p); %% rad to deg conversion
f1 = ((0:length(Y)-1)*Fs/length(Y))\';
figure
subplot(2,1,1)
plot(f1,m)
title(\'Magnitude\')

subplot(2,1,2)
plot(f1,p)
title(\'Phase\')

%%% saving text file in ASCII foramt%%
% 1st column- frequency%%
% 2nd column- magnitude components%%
% 3rd column- phase (deg) components%%
A(:,1)=f1;
A(:,2)=m;
A(:,3)=p;
save(\'fre_amp_phase.txt\',\'A\', \'-ASCII\')