Matlab Fourier

fourier.m

%% Matlab-Skript zur Vorlesung Informationsuebertragung am 14.11.2014
% Autor: Valentin Burger, Michael Seufert, Steffen Gebert

clear all; %set(0, 'defaultlinelinewidth', 1);
figure(1);clf;box on;
A=5;                        % Amplitude
w=4;                        % Frequenz
phi=pi/2;                   % Phasenverschiebung
x=(0:0.001:2)*pi;
y=A*sin(x*w+phi);
plot(x,y);
set(gca,'XTick',(0:0.5:2)*pi)
set(gca,'XTickLabel',{'0' 'pi/2' 'pi' '3 pi/2' '2 pi'})
xlabel('$x$','Interpreter','latex');
ylabel('$A\cdot sin(x\cdot\omega + \phi)$','Interpreter','latex');
xlim([0 2*pi])

figure(2);clf;box on;
f=4;                        % Frequenz (der Schwingung) in Hertz
phi = 0;                    % Phasenverschiebung
Fs = 1000;                  % Abtastfrequenz
t = 0:1/Fs:1;               % Abtastzeitpunkte
a = A*sin(t*2*pi*f+phi);
plot(t,a);
xlabel('Zeit t / s');ylabel('Amplitude');
%% (2) Ueberlagerung zweier Schwinungen
figure(1);clf;hold all;box on;
Fs = 1000;
t = 0:1/Fs:1;
A=[1 0.5];
f=[1 2];
phi=[0 0];
plot(t,A(1)*sin(t*f(1)*2*pi+phi(1)),'red', 'LineWidth', 2);
plot(t,A(2)*sin(t*f(2)*2*pi+phi(2)),'black', 'LineWidth', 2);
y=A(1)*sin(t*f(1)*2*pi+phi(1));
y=y+A(2)*sin(t*f(2)*2*pi+phi(2));
plot(t,y,'blue', 'LineWidth', 2)
xlabel('Zeit t / s');ylabel('Amplitude');
legend('signal1','signal2','signal(1+2)')
%% (3) Uberlagerung mehrerer Schwingungen
A=[2 1 1 0.5];
f=[2 4 7 10];
phi=[0 0 pi -pi/2];
Fs = 1000;
t = 0:1/Fs:1;
%Plot der einzelnen Schwingungen
%figure(1);clf;box on;hold all;
%for i=1:4 plot(t,A(i)*sin(t*f(i)*2*pi+phi(i)));end
figure(2);clf;box on;
y=A(1)*sin(t*f(1)*2*pi+phi(1));
y=y+A(2)*sin(t*f(2)*2*pi+phi(1));
y=y+A(3)*sin(t*f(3)*2*pi+phi(1));
y=y+A(4)*sin(t*f(4)*2*pi+phi(1));
plot(t,y)
xlabel('Zeit t / s');ylabel('Amplitude');
%% (4) FFT komplettes Spektrum
figure(1);clf;box on;
f=4;
A=5;
phi=0;
Fs = 1000;
t = 0:1/Fs:1;
L = length(t);
a = A*sin(t*2*pi*f+phi); 
plot(t,a);
xlabel('Zeit t / s');ylabel('Amplitude');

figure(2);clf;box on;
% matlab fft-Funktion shiftet Spektrum und multipliziert mit L
ft = fftshift(fft(a)/L);
% man erhaelt L Werte im Frequenzspektrum, die von -Fs/2 bis Fs/2 reichen
freq = Fs/2*linspace(-1,1,L);
% hier wird der Betrag geplotet, vgl Funktionen fuer Real/Imaginaerteil
% real/imag
plot(freq, abs(ft),'o','MarkerSize',10)
xlim([-10 10])
xlabel('Frequenz f / Hz');ylabel('Amplitude');
%% (5) Darstellung des Spektrums in positivem Frequenzbereich
% diese Zelle wird in fourierS-Funktion zusammengefasst
figure(1);clf;box on;
f=4;
A=5;
phi=0;
Fs = 1000;
t = 0:1/Fs:1;
L = length(t);
a = A*sin(t*2*pi*f+phi);
% matlab fft-Funktion shiftet Spektrum und multipliziert mit L
FFTa = fft(a)/L;
% wir nehmen L/2 bzw. (L-1)/2 Werte im positivem Spektrum von 0 bis Fs/2
freq = Fs/2*linspace(0,1,(L-1)/2);
% Amplitude setzt sich aus positivem und negativem Frequenzbereich zusammen,
% also mal 2
plot(freq, 2*abs(FFTa(1:(L-1)/2)),'o','MarkerSize',10);
xlabel('Frequenz in Hz / s');ylabel('Amplitude');
xlim([0 10])
% %% (6) Kopie von oben
% figure(1);clf;box on;hold all;
% xlabel('Zeit t / s');ylabel('Amplitude');
% figure(2);clf;box on;hold all;
% xlabel('Frequenz f / Hz');ylabel('Amplitude');xlim([0 10])
% Fs = 1000;
% T = 1/Fs;
% t = 0:T:1;
% A=[1 0.5];
% f=[1 2];
% phi=[0 0];
% a1 = A(1)*sin(t*f(1)*2*pi+phi(1));
% figure(1); plot(t,a1,'--r');
% a2 = A(2)*sin(t*f(2)*2*pi+phi(2));
% figure(1); plot(t,a2,'--black');
% a12=a1 + a2;
% figure(1); plot(t,a12,'blue')
% legend('signal1','signal2','signal(1+2)')
% %% (7) In Ueberlagerten Schwingung enthaltene Frequenzen
% [FTa1, freq] = fourierS(a1,Fs);
% figure(2); plot(freq,abs(FTa1),'xr','MarkerSize',10)
% 
% [FTa2, freq] = fourierS(a2,Fs);
% figure(2); plot(freq,abs(FTa2),'oblack','MarkerSize',10)
% 
% a12=a1 + a2;
% 
% figure(1); plot(t,a12,'blue')
% 
% [FTa12, freq] = fourierS(a12,Fs);
% figure(2); plot(freq,abs(FTa12),'+b','MarkerSize',10)
% 
% figure(1);legend('signal1','signal2','signal(1+2)')
% figure(2);legend('signal1','signal2','signal(1+2)')
% %% (8) bisher betrachten wir nur den Betrag, wir unterschlagen die Phase
% figure(1);clf;box on;hold all;
% xlabel('Zeit t / s');ylabel('Amplitude');
% figure(2);clf;box on;hold all;
% xlabel('Frequenz f / Hz');ylabel('Amplitude');xlim([0 10]);title('Betrag')
% figure(3);clf;box on;hold all;
% xlabel('Frequenz f / Hz');ylabel('Amplitude');xlim([0 10]);title('Realteil')
% figure(4);clf;box on;hold all;
% xlabel('Frequenz f / Hz');ylabel('Amplitude');xlim([0 10]);title('Imaginaerteil')
% Fs = 1000;
% t = 0:1/Fs:1;
% A=[1 1];
% f=[3 3];
% phi=[0 pi/2];
% a1 = A(1)*cos(t*f(1)*2*pi+phi(1));
% figure(1); plot(t,a1,'--r');
% [FTa1, freq] = fourierS(a1,Fs);
% figure(2); plot(freq,abs(FTa1),'or','MarkerSize',10)
% figure(3); plot(freq,real(FTa1),'or','MarkerSize',10)
% figure(4); plot(freq,imag(FTa1),'or','MarkerSize',10)
% 
% a2 = A(2)*cos(t*f(2)*2*pi+phi(2));
% figure(1); plot(t,a2,'--black');
% [FTa2, freq] = fourierS(a2,Fs);
% figure(2); plot(freq,abs(FTa2),'+black','MarkerSize',10)
% figure(3); plot(freq,real(FTa2),'+black','MarkerSize',10)
% figure(4); plot(freq,imag(FTa2),'+black','MarkerSize',10)
%% (9) Rechtecksimpuls
Fs = 100;
t = -1:1/Fs:1;
a=zeros(1,length(t));
a(t>-0.5 & t<0.5) = 1;

figure(1);clf;box on;
plot(t,a, 'LineWidth', 2)
ylim([0 1.5])

%% (10) animation

j=2;
x = [zeros(1, 100), j*ones(1, 100), zeros(1, 100)];
y = fftshift(fft(x, 1000));

f1=figure(1);
hold on;

for i = 1: 100
    clf
   xi = abs(ifft(y(length(y)/2-((i-1)*10)/2:length(y)/2+((i-1)*10)/2), 1000));

   plot(xi(1:300), 'LineWidth', 2);    

   xlabel('Zeit t / s');
   ylim([0, j+0.5]);
   refresh(f1);
 pause(0.1);
   M(i) = getframe;
   %clf;

end

 %plot(xi(1:300), 'LineWidth', 2);   ylim([0, 1.1]);
 


% koeffizienten
[FTa, freq, FFTa] = fourierS(a,Fs);
figure(2);clf;box on;
plot(freq,real(FTa),'o','LineWidth',2,'MarkerEdgeColor','k',...
                'MarkerFaceColor','r','MarkerSize',10)
xlim([0 50]);ylim([-1 1]);
xlabel('f/Hz');ylabel('Amplitude');


%% WAV-Datei einlesen und Frequenzspektren plotten
figure(3);clf;

[a, Fs] = wavread('run.wav');
t = (0:length(a)-1)'/Fs;

[FTa, freq, FFTa] = fourierS(a,Fs);

subplot(2,2,1), plot(t,a), xlabel('t/s'), ylabel('Amplitude'), title('run.wav');
subplot(2,2,2), plot(freq, abs(FTa)), xlabel('f/Hz'), ylabel('Amplitude');

[a, Fs] = wavread('ring.wav');
t = (0:length(a)-1)'/Fs;

[FTa, freq, FFTa] = fourierS(a,Fs);

subplot(2,2,3), plot(t,a), xlabel('t/s'), ylabel('Amplitude'), title('ring.wav');
subplot(2,2,4), plot(freq, abs(FTa)), xlabel('f/Hz'), ylabel('Amplitude');
%% Tiefpass

[a, Fs] = wavread('2.wav');
figure(1);clf;box on;hold all;
xlabel('Zeit t / s');ylabel('Amplitude');
t = (0:length(a)-1)'/Fs;
plot(t,a)

[FTa, freq, FFTa] = fourierS(a,Fs);
figure(2);clf;box on;hold all;
xlabel('Frequenz f / Hz');ylabel('Amplitude');xlim([0 10]);

plot(freq,abs(FTa),'blue')
xlim([0 15000]);
%%
figure(2);
plot(freq,0.1.*ones(1,length(freq)).*(freq<=4000),'-g')
%%
% geshifteter Rechtecksimpuls von -4000 Hz bis 4000 Hz
index = find(freq>4000,1);
rect4000 = ones(1,length(FFTa))';
rect4000(index:end-index) = 0;
% Tiefpassfilter (in FFTa is geshiftetes Spektrum)
cleanFFTa = FFTa.*rect4000;

% Ruecktransformation des gefilterten Spektrums
cleana = ifft(cleanFFTa);
[cFTa, freq] = fourierS(cleana,Fs);
figure(2);plot(freq, abs(cFTa),'r')

figure(1);plot(t,cleana,'r')
% WAV-Datei schreiben
wavwrite(cleana,Fs,'2clean.wav')
%% subplots
figure(3);clf;
subplot(2,2,1), plot(t,a), xlabel('t/s'), ylabel('Amplitude'), title('gestoertes Signal');
subplot(2,2,2), plot(freq, abs(FTa)), xlim([0 8100]), ylim([0 0.2]),...
    hold on, plot(freq,0.1*ones(1,length(freq))*(freq<=4000),'-r'),...
    xlabel('f/Hz'), ylabel('Amplitude');
subplot(2,2,3), plot(t,cleana), xlabel('t/s'), ylabel('Amplitude'), title('Signal nach Tiefpass');
subplot(2,2,4), plot(freq, abs(cFTa)), xlim([0 8100]), xlabel('f/Hz'), ylabel('Amplitude');


%% Ende