Leistedichtespektrum eines abgetasteten Signals

[Ankaios A. (ankaios)]

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!

    clear all;
    %% Configuration
    T               = 2^10;     % [s]
    f_max           = 2;        % [Hz]
    fs              = f_max*4;  % [Hz] Nyquist–Shannon sampling theorem: fs>=f_max*2
    SampleLength    = 16;        % [-]
    %% Analytical Spectra 
    df              = 1/T;
    f               = [df:df:f_max];
    n_FFT           = fs/df;
    dt              = 1/fs;
    t               = 0:dt:T-dt;
    %-----------
    S               = 1./(1+10*f);
    %-----------
    RECT           = sinc(f*SampleLength*dt);
    S_MAV           = S .*(RECT.*conj(RECT));
    %% Time Signals
    phi             = rand(size(f))*2*pi;   % random phase shift
    A               = sqrt(S/(1/2/df));     % Amplitudes
    % time signal is a sum of all frequencies
    x               = (A*cos(2*pi*f'*t+repmat(phi',1,n_FFT)));
    % Moving AVerage
    x_MAV          = filter(ones(1,SampleLength)/SampleLength,1,x);
    % Zero Order Hold
    x_ZOH            = reshape(repmat(x(1:SampleLength:end),SampleLength,1),1,[]);
    %% Spectra from data
    nBlocks             = 16;
    nDataPerBlock        = floor(n_FFT/nBlocks/2)*2; % should be even
    vWindow            = hann(nDataPerBlock);
    %-----------
    [S_D        f_D]   = pwelch(detrend(x,'constant'),     vWindow,[],[],fs);
    [S_D_MAV    f_D]    = pwelch(detrend(x_MAV,'constant'), vWindow,[],[],fs);
    [S_D_ZOH    f_D]    = pwelch(detrend(x_ZOH,'constant'), vWindow,[],[],fs);   
    %% Plots
    figure('Name', 'Time')
    hold on;grid on;box on;
    plot(t,x,'.-');
    plot(t,x_MAV ,'g');
    plot(t,x_ZOH ,'o-r');
    xlim([0 8])
    xlabel('t [s]')
    ylabel('x [-]')
    legend('x','x_{MAV}','x_{ZOH}')
    %-----------
    figure('Name', 'Spectra')
    hold on;grid on;box on;
    plot(f_D,S_D,'LineWidth',3)
    plot(f_D,S_D_MAV,'g','LineWidth',3)
    plot(f_D,S_D_ZOH,'r','LineWidth',3)
    plot(f,S,'c','LineWidth',1)
    plot(f,S_MAV,'k','LineWidth',1)
    legend('S_{D}','S_{D,MAV}','S_{D,ZOH}','S_{A}','S_{A,MAV}')
    set(gca,'XScale','log')
    set(gca,'YScale','log')
    ylabel('S [1/Hz]')
    xlabel('f [Hz]')
    xlim([df*nBlocks f_max])
    ylim([1e-5 1e1])

[Herr Shannon (Gast)]

Ankaios Argo schrieb:
> Bisher sind alle meine
> Versuche gescheitert, hier eine "analytische" Lösung basierend auf der
> Gleichung des original Spektrums zu finden.

Tippe auf schwer unterabgetastet.

[Aasgeier Sigurvinson (Gast)]

Ankaios Argo schrieb:
> ich einfach alle Frequenzen
>
> überlagert habe.

das wäre dann aber gleichspannung

kann nicht sein

[Ankaios A. (ankaios)]

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:

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.

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:

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

also

S_ZOH           = S .*???;