Leistedichtespektrum eines abgetasteten Signals

Ankaios A. OP #3169498
Lesenswert?
• ▲
▼ 
Guten Abend!

Gerade versuche ich, das Spektrum eines Signals zu verstehen, welches 
mit einem "zero-order-hold" abgetastet wird.

Dazu habe ich ein Zeitsignal (code und plots siehe unten) aus einem 
analytischen Spektrum gebaut, indem ich einfach alle Frequenzen 
überlagert habe. Aus dem Zeitsignal konnte ich dann recht gut das 
original Spektrum wieder rekonstruieren.
Danach habe ich das Zeitsignal mit einem laufenden Mittelwert gefiltert 
und auch hier hat das Spektrum aus den Daten mit der Theorie recht gut 
übereingestimmt: Da die Faltung im Zeitbereich zu einer Multiplikation 
im Frequenzbereich wird, musste ich nur das ursprüngliche Spektum mit 
sinc*conj(sinc) multiplizieren.

Danach habe ich aus dem original Zeitsignal jeden 16ten Wert genommen 
und gehalten und auch das Spektrum begrechnet. Bisher sind alle meine 
Versuche gescheitert, hier eine "analytische" Lösung basierend auf der 
Gleichung des original Spektrums zu finden.

Hat da jemand eine Idee?

Nach Folie 9 von
http://www.cis.hut.fi/Studies/T-61.140/Luennot/Ch7.pdf

Sollte die Fourier-Transformierte eine unendliche Summe von verschobenen 
Fourier-Transformierten sein. Aber wie setzt man das konkret um?

Mit für jeden Hinweis/Idee sehr dankbar!

1
    clear all;

2
    %% Configuration

3
    T               = 2^10;     % [s]

4
    f_max           = 2;        % [Hz]

5
    fs              = f_max*4;  % [Hz] Nyquist–Shannon sampling theorem: fs>=f_max*2

6
    SampleLength    = 16;        % [-]

7
    %% Analytical Spectra 

8
    df              = 1/T;

9
    f               = [df:df:f_max];

10
    n_FFT           = fs/df;

11
    dt              = 1/fs;

12
    t               = 0:dt:T-dt;

13
    %-----------

14
    S               = 1./(1+10*f);

15
    %-----------

16
    RECT           = sinc(f*SampleLength*dt);

17
    S_MAV           = S .*(RECT.*conj(RECT));

18
    %% Time Signals

19
    phi             = rand(size(f))*2*pi;   % random phase shift

20
    A               = sqrt(S/(1/2/df));     % Amplitudes

21
    % time signal is a sum of all frequencies

22
    x               = (A*cos(2*pi*f'*t+repmat(phi',1,n_FFT)));

23
    % Moving AVerage

24
    x_MAV          = filter(ones(1,SampleLength)/SampleLength,1,x);

25
    % Zero Order Hold

26
    x_ZOH            = reshape(repmat(x(1:SampleLength:end),SampleLength,1),1,[]);

27
    %% Spectra from data

28
    nBlocks             = 16;

29
    nDataPerBlock        = floor(n_FFT/nBlocks/2)*2; % should be even

30
    vWindow            = hann(nDataPerBlock);

31
    %-----------

32
    [S_D        f_D]   = pwelch(detrend(x,'constant'),     vWindow,[],[],fs);

33
    [S_D_MAV    f_D]    = pwelch(detrend(x_MAV,'constant'), vWindow,[],[],fs);

34
    [S_D_ZOH    f_D]    = pwelch(detrend(x_ZOH,'constant'), vWindow,[],[],fs);   

35
    %% Plots

36
    figure('Name', 'Time')

37
    hold on;grid on;box on;

38
    plot(t,x,'.-');

39
    plot(t,x_MAV ,'g');

40
    plot(t,x_ZOH ,'o-r');

41
    xlim([0 8])

42
    xlabel('t [s]')

43
    ylabel('x [-]')

44
    legend('x','x_{MAV}','x_{ZOH}')

45
    %-----------

46
    figure('Name', 'Spectra')

47
    hold on;grid on;box on;

48
    plot(f_D,S_D,'LineWidth',3)

49
    plot(f_D,S_D_MAV,'g','LineWidth',3)

50
    plot(f_D,S_D_ZOH,'r','LineWidth',3)

51
    plot(f,S,'c','LineWidth',1)

52
    plot(f,S_MAV,'k','LineWidth',1)

53
    legend('S_{D}','S_{D,MAV}','S_{D,ZOH}','S_{A}','S_{A,MAV}')

54
    set(gca,'XScale','log')

55
    set(gca,'YScale','log')

56
    ylabel('S [1/Hz]')

57
    xlabel('f [Hz]')

58
    xlim([df*nBlocks f_max])

59
    ylim([1e-5 1e1])
Angehängte Dateien:
Antwort Zitieren Bearbeiten Thread verschieben Thread sperren Anmeldepflicht aktivieren Anpinnen Thread löschen Thread mit anderem zusammenführen Beitrag melden
Ankaios A. OP #3169581
Lesenswert?
• ▲
▼ 
Vielen Dank für die Antworten!

Ich bin ein rechter Neuling in der Signalverarbeitung, aber ich dachte, 
dass ich das Signal schon richtig abgetastet habe mit einer Frequenz, 
die viermal so hoch ist wie die höchste Frequenz im Signal:

1
fs              = f_max*4;


Vielleicht hab ich mich falsch ausgedrückt, aber das mit den 
überlagerten Frequenzen war so gemeint, dass ich aus dem Spektrum die 
Amplituden genommen habe und dann ein Zeitsignal zu jeder Frequenz mit 
der zugehörigen Amplitude mit einem zufälligen Phasenversatz erzeugt 
habe und dann einfach addiert. Wenn ich daraus dann das 
Leistungsspektrum berechne, passt es wieder zu meinem ursprünglichen 
Spektrum.

1
x               = (A*cos(2*pi*f'*t+repmat(phi',1,n_FFT)));


Was ich gerne hätte, wäre eine analytische Lösung wie die zu dem 
gleitenden Mittelwert:

1
S_MAV           = S .*(sinc(f*T_MAV).*conj(sinc(f*T_MAV)));


also

1
S_ZOH           = S .*???;
Antwort Zitieren Bearbeiten Löschen Beitrag melden
Thread beobachten

Antwort schreiben

Bitte melde dich an, um einen Beitrag zu schreiben.

Bestehender Account

oder

Mit Google-Account einloggen

Noch kein Account?

Die Registrierung ist kostenlos und dauert nur eine Minute.

Jetzt registrieren