How to visualize this statement regarding Conjugate Symmetry A property of real signals states that if $x(t)$ is real then the Fourier series coefficient (frequency spectrum) is given by :
$$
c_{k}=c_{-k}^{*}
$$
In polar form, this can be expressed as :
$$
|c_{k}|=|c_{-k}|\qquad\text{and}\qquad\arg c_{k}=-\arg c_{-k}
$$
It appears here that the magnitude spectrum obeys even symmetry and the phase spectrum obeys odd symmetry. It is said in my textbook :

Due to symmetry in the frequency spectra of real signals, we typically ignore negative frequencies when dealing with such signals.
In the case of signals that are complex but not real, frequency spectra do not possess the above symmetry,
and negative frequencies become important.

For $(1)$, how can I visualize this to see why the negative frequency components do not matter anymore$?$, does this mean the Fourier series representation of $x(t)$ becomes :
$$
x(t)=\sum_{k=0}^{\infty}c_{k}e^{jk\omega_{0} t}\;???
$$
For $(2)$, is this statement true because the complex exponentials in the frequency domain will scale a given signal to a decaying or growing signal?

For $(1)$, how can I visualize this to see why the negative frequency components do not matter anymore$?$, does this mean the Fourier series representation of $x(t)$ becomes :
$$
x(t)=\sum_{k=0}^{\infty}c_{k}e^{jk\omega_{0} t}\;???
$$

No, the Fourier series is still what it was before this conjugate symmetry nonsense reared its ugly head. But what happens is that
\begin{align}
x(t) &= \sum_{k=-\infty}^{\infty}c_{k}e^{jk\omega_{0} t} & {\scriptstyle\text{original definition}}\\
&= \sum_{k=-\infty}^{-1}c_{k}e^{jk\omega_{0} t} + c_0 + \sum_{k=1}^{\infty}c_{k}e^{jk\omega_{0} t}& {\scriptstyle\text{breaking up the sum}}\\
&= \sum_{k=1}^{\infty}c_{-k}e^{-jk\omega_{0} t} + c_0 + \sum_{k=1}^{\infty}c_{k}e^{jk\omega_{0} t}& {\scriptstyle\text{changing indices}}\\
&= c_0 + \sum_{k=1}^{\infty}c_{k}e^{jk\omega_{0} t} + c_{-k}e^{-jk\omega_{0} t} & {\scriptstyle\text{combining sums}}\\
&= c_0 + \sum_{k=1}^{\infty}c_{k}e^{jk\omega_{0} t} + c_{k}^*e^{-jk\omega_{0} t} & {\scriptstyle c_k = c_{-k}^* \implies c_{-k} = c_k^*}\\
&= c_0 + \sum_{k=1}^{\infty}c_{k}e^{jk\omega_{0} t} + \left(c_{k}e^{jk\omega_{0} t}\right)^* & {\scriptstyle a^*b^* = (ab)^*}\\
&= c_0 + 2\sum_{k=1}^{\infty}\operatorname{Re}\left(c_{k}e^{jk\omega_{0} t}\right)
\end{align}
showing that one needs to compute or know the $c_k$ values for nonnegative $k$ only; the values of $c_k$ for $k < 0$ are not needed. Anal-retentives who believe in crossing their eyes and dotting their teas can always determine the values of the $c_{-k}$ by taking the complex conjugates of the known $c_k$'s.


Is it possible to specify number of taps and cutoff frequency independently? When specifying the properties of a moving average filter, is it possible to specify the number of taps and the cutoff frequency?
My guess is no since taps (the number of samples in a moving average filter) is what determines the cutoff frequency according to the frequency response equation.
If I'm wrong, how is it possible to specify the two?

If you talk about a moving average filter then you only have a single degree of freedom, namely the filter length (i.e., the number of taps). Consequently, given the number of taps, the magnitude of the frequency response is fixed and the cut-off frequency is given.
Note that you could also define the delay of the moving average filter as another degree of freedom, but that wouldn't change the filter's magnitude response and its cut-off frequency.
For general FIR filters with $N$ taps you have $N$ degrees of freedom, so you can clearly choose the cut-off frequency independently of the number of taps. But this is not the case for a moving average filter.


How to form Kalman filtering matrices for a problem with variable acceleration? Assuming we have

time vector $T$ with constant time step $dt$
position vector $X$
velocity vector $V$
acceleration vector $A$

All vectors $X, V, A$ have noise on their measurement ( $n_x$ , $n_v$ , $n_a$ ). And all should be filtered.

How to form Kalman filter for the the matrices:

$$ x_k= A x_{k-1} + B u_k + w_{k-1} $$
$$ y_k= H x_k + v_k $$

where $w$ and $n$ are normal noise. Each $x_k$ is $\begin{bmatrix} X_k  & V_k & A_k\end{bmatrix}^T$

I came up with two matrix assumption for kalman filtering:

which one is correct?
$$A=\begin{bmatrix} 1 & dt & \frac{1}{2}dt^2 \\ 0 & 1 & dt \\ 0 & 0 & 1 \end{bmatrix}$$
$$B=\begin{bmatrix}0 & 0 & 0\end{bmatrix}^T$$
$$H=I$$

or this one?

$$A=\begin{bmatrix} 1 & dt & \frac{1}{2}dt^2 \\ 0 & 1 & dt \\ 0 & 0 & 0 \end{bmatrix}$$
$$B=\begin{bmatrix} 0 &  0 & 1 \end{bmatrix}^T$$
$$H=I$$

The acceleration is not constant and there is a variable input acceleration. I need to fix acceleration as well as velocity and position.

Please explain why which one is correct and the other is wrong.

This is my code. How to extend it for variable acceleration?

time_resolution=0.016;
tmp1=diag((1:9)*0+1,3);
tmp2=diag((1:9)*0+1,6);
A_T=eye(9)+tmp1(1:9,1:9)*time_resolution+0.5*tmp2(1:9,1:9)*time_resolution*time_resolution;
A=A_T;
C=eye(size(A));
p_tmp=[position_x;position_y;position_z];
v_tmp=[velocity_x;velocity_y;velocity_z];
a_tmp=[acceleration_x;acceleration_y;acceleration_z];
X=[p_tmp;v_tmp;a_tmp];

xh_=X*0;
yh_=X*0;
inov=X*0;
xh=zeros(size(X)+[0,1]);
xh(:,1)=X(:,1);

R=1;
Q=100*diag([.1 .1 .1 .5 .5 .5 10 10 10]);
Px=eye(size(A));

for i = 1 : size(X,2)
    % x(t|t-1) = A*x(t-1|t-1)
    xh_(:,i) = A * xh(:,i);

    % P(t|t-1) = A*P(t-1|t-1)*A' + Q
    Px_ = A*Px*A' + Q;

    % K(t) = P(t|t-1) * C' * inv(C*P(t|t-1) * C' + R) )
    K = Px_ * C' /(C*Px_*C' + R);

    % y(t|t-1) = C*x(t|t-1) + R
    yh_(:,i) = C * xh_(:,i) + R;

    % innovation error  y(t) - y(t|t-1)
    inov(:,i) = X(:,i) - yh_(:,i);

    % x(t|t) = x(t|t-1) + K(t)*(y(t)-y(t|t-1))
    xh(:,i+1) = xh_(:,i) + K * inov(:,i);

    % P(t|t) = (I - K(t)*C) * P(t|t-1)
    Px = Px_ - K*C*Px_;
end

Take the first model and train your model noise matrix $\bf{Q}$ to account for the variable acceleration. A popular model is the Singer model from 1990 (see [1]).

[1] Singer, R. A., “Estimating Optimal Tracking Filter Performance for Manned Maneuvering Targets,” IEEE Transactions on Aerospace and Electronic Systems, vol. AES-6, no. 4, .pp. 473-483, July 1990.


Pole-Zero Plots Diagram Explanation I experiment with the coefficient quantization in an IIR filter. I change some values like the passband ripple and also how many bits i want for the quantization method. In the first picture i have 1db passband ripple , 'round'and 7-bits as the quantization method:

In the second picture i have 1.3db as the passband ripple ,'round' and 6-bits as quantization method.

These are the PZP diagrams. I dont understand how i can extract useful information from the diagrams in order to say which design is better or to make comparisons between these 2 designs.How the Poles and Zeros affect my decision for what is better?

edit: the images becomes bigger and more visible when you copy the url and open or select view image. Also with the green is the quantized and with blue the un-quantized coefficients.

edit2:As Matt asks for magnitude frequency response I have add the extra information:

1db passband ripple , 'round'and 7-bits

1.3db as the passband ripple ,'round' and 6-bits

EDIT: After the answer from Matt i need to add more information:

I want to find a set of quantized coefficients that meet the specification with the minimum number of bits.
 My specification is an 8th-order IIR filter with the following transfer function magnitude response:
\begin{align*}
 0\,\mathrm{db} \pm 1.5 \,\mathrm{db} & \quad \text{for } \quad 0.2<|v|<0.3 \cr
 \lt −60\,\mathrm{db} & \quad\text{for } \quad |v|<0.14\text{ and } |v|>0.36
\end{align*}

I use this matlab function:  $ [b,a]=ellip(n,Rp,Rs,Wp) , $

after some trials i found that (the filter is implemented using second-order sections) :  $ [b,a]=ellip(4,1,60,2*[0.2,0.3]) $  and using 7-bit quantization i take the diagram above that i think verify the specification (if i make a mistake please correct me) ,

but if a use: $ [b,a]=ellip(4,1.3,60,2*[0.2,0.3]) $ and 6-bit quantization also i take a reasonable result. My concern is if my thought and procedure is correct or i have a mistake and i need to stick with 7-bits quantization?

It is indeed hard to compare the two designs just using the pole-zero diagrams. I think it would be much more useful to have a look at the (magnitude of the) frequency responses, and compare the change of the frequency response before and after coefficient quantization. The only thing I can say from looking at the pole-zero plots is that you have designed an eighth order bandpass filter with center frequency $f_s/4$ (by $f_s$ I mean the sampling frequency).


Does dividing the magnitude spectrum of white noise by sqrt(2) give an RMS magnitude spectrum? I understand that the RMS Amplitude of a sinusoidal signal is around 0.707 ($\frac1{\sqrt2}$) times the Peak Value, but this is not true for noise.

However, an FFT of a noise signal indicates magnitude in bins at equally spaced frequencies, which means the magnitude represented in each bin, is the magnitude of a periodic signal (sinusoid).

So, I assume that multiplying the magnitude of an FFT by 0.707, should give me RMS magnitude spectrum on the Y-axis, instead of the peak value. Am I correct in assuming this, or am I missing something?

However, an FFT of a noise signal indicates magnitude in bins at equally spaced frequencies, which means the magnitude represented in each bin, is the magnitude of a periodic signal (sinusoid).

No; it just means that when projecting the noise signal onto a sinusoidal signal, you don't get nothing.

Assuming your noise realization in each is independent identically distributed (usually the way we model noise!), then the (rotated) sum of such noise samples approximates a normal distribution, no matter how the (finite-variance) original distribution looked like. The magnitude in each bin is an approximation (because by definition, your DFT has finite length!) to the noise power density; in an uncorrelated noise process, that's the variance==power==rms².


why does butterworth IIR filter for a square pulse show ripple on edges in time domain but not the same for sine wave I am using butterworth filter for both square pulse and sine wave. In case of square pulse, butterworth produces some sharp ripple of edges of square pulse in time domain for varying butterworth order (wn > 2). But in case of sine wave, I don't see any ripple. Could anyone please explain the reason behind that? Matlab code is following along with plots.
close all

clear all
T = 1/10;
t = linspace(0,T,1001);

s = square (2*pi*100*t);

fs=10/200; % sampling frequency
fc=0.2; % fc normalized
offset=0;
amp=3.3;
duty=50;
t=0:0.01:100;%100 seconds
sq_wav=offset+amp*square(2*pi*fs.*t,duty);
sin_wav = offset + amp*sin(2*pi*fs.*t);
disp(length(sin_wav));

N = length(sq_wav);
disp(N);
#disp(t(1:10));
t1 = [0:N-1];
disp(t1(1:10));
disp(fs);
wn = 15; % filter order

[b,a] = butter(wn, fc);

[h,w] = freqz(b,a);

out_sq = filter(b,a,sq_wav);
out_sin = filter(b,a,sin_wav);

subplot (4, 1, 1)
plot(t(1:end),sq_wav(1:end))
title("square pulse");

subplot (4, 1, 2)
plot(t(1:end),sin_wav(1:end))
title("sine pulse");

subplot (4, 1, 3)
plot(t(1:end),out_sq(1:end))
title("square pluse after butterworth");

subplot (4, 1, 4)
plot(t(1:end),out_sin(1:end))
title("sine pluse after butterworth");

Here is the output of this code:

This happens because when lowpass filtering, you remove the high-frequency content that gives the square pulses their sharp edges. Now, if you're not yet versed in Fourier analysis of signals this will be a bit tough, but I'll provide some graphs to show what I mean.
Let's take a simple pulse and a sinusoid like the ones you have. We're going to filter them with a simple 4-th order Butterworth filter. Below are the time-domain and frequency-domain plots of the signals. The spectrum of the filter is also included. The spectrum plots are in decibels (dB).

If you don't know, filtering in the time domain is multiplication in the frequency domain. So you can look at the frequency spectrums and multiply them together, seeing where the attenuation takes place. You can see that the pulse has more significant contributions from high frequencies compared to the sinusoid. So after filtering, the sinusoid will be less distorted than the pulse which is exactly what we see below

We can take this a step further and remove more frequency content by using a tighter filter. Below are the same spectrums with the new filter. You see that the pulse is now even more distorted, while the sinusoid remains relatively unchanged.


XAudio2 filter radian frequency In the documentation for the [XAudio2 filter parameters] it mentions

Filter radian frequency calculated as (2*sin(pi*(desired filter cutoff frequency)/sampleRate))

Where does this formual come from? What is "Filter radian frequency" and "filter cutoff frequency" and what is the difference. Can someone translate this terminology to mathematics of DFT?

That is simply converting analog frequency into digital frequency. For example, if you want a cutoff of 1kHz at a sample rate of 44.1kHz, the corresponding digital frequency is 2*pi*1000/44100 = 0.142... radiants.
In terms of DFT the digital frequency is the $ \omega $ in the $e^{^-j\cdot \omega }$ term of the DFT


What is the cut-off frequency for this low pass filter? I don't have any background in analog lowpass filters; and I am curious how can i decode this circuit.

Look for Sallen and Key filter.
Fc = 1/(2 * pi * sqrt(R1 * R2 * C1 * C2))
Or use a calculator, found one here: Sallen-Key Low-pass Filter Design Tool
(It appears to handle "k" and "n", as a bonus)
Fc = 232 Hz, if you want to check your math.


Signal power without FFT I have an EEG signal thats about 1second long, sampling frequency of 256Hz, measured in microvolts, and I want to calculate the power of a specific frequency band (8-12Hz)

I believe the correct way to do this is to do an FFT which will give you frequency domain information, from which you can directly extract power. Is that correct?

I'm looking for a way to do this without FFT.

I have applied a butterworth filter to my signal, bandpass at 8-12Hz. This, I think, will result in a filtered signal that shows only the contribution from sinusoids within those frequency ranges.

Is there a way to extract power directly from this bandpassed signal without having to run an FFT? Would it simply be the average (squared?) voltage across the entire signal? Or would it be more complex than that?

I'm doing this in Python, so if there's any specific code I can take a look at then that'd be useful

After performing the band-pass filter, the energy in that frequency band is given by:

sum( abs( signal ) ^ 2 )

By abs I mean the magnitude of a complex number, or the absolute (positive) value of a real number.

This works because of Plancherel's theorem.


Mathematical equation for the sound wave that a piano makes Is there a mathematical function that can output the sound that a piano makes at different frequencies?

So if the simplest equation for a sound wave (at a given frequency and for a given sample rate) is

$ sin(\frac{2 \cdot \pi \cdot x \cdot f}{samplerate}) $

where f is frequency, is there a similar equation for the sound that a piano makes?

EDIT: I'm not looking for anything too complicated here, just a wave that approximates this estimation.   I've been playing around on desmos to see if I can get something similar:

This is what Guest's wave sounds like... not much like a piano.

Piano tones are actually quite complicated.  The initial strike is full of non-harmonic tones that quickly dampen down.  After that, the harmonics vary from note to note, piano to piano.  I suggest you look at some actual sample files, capture a whole number of somewhat consistent waveforms, do a DFT, and read the coefficients from there.

As for the other answer that was given, some one may find this interesting:

$$ \sin^3\left(\pi x\right) + \sin \left( \pi \left( x + \frac{2}{3} \right) \right) $$

$$ = \left( \frac{ e^{i\pi x} - e^{-i\pi x} }{2i}  \right)^3 + \sin\left( \pi x+ \frac{2}{3}\pi \right) $$

$$ = \frac{ e^{i3\pi x} - 3 e^{i\pi x} + 3 e^{-i\pi x} - e^{-i3\pi x} }{-8i} + \sin( \pi x ) \cos \left(\frac{2}{3}\pi \right) + \cos( \pi x ) \sin \left(\frac{2}{3}\pi \right) $$

$$ = -\frac{ 1 }{4} \sin(3\pi x) + \frac{ 3 }{4} \sin(\pi x)  + \sin( \pi x ) \left(\frac{-1}{2} \right) + \cos( \pi x ) \left(\frac{\sqrt{3}}{2} \right) $$

$$ = -\frac{ 1 }{4} \sin(3\pi x) + \frac{ 1 }{4} \sin(\pi x)   + \frac{\sqrt{3}}{2} \cos( \pi x )  $$

Hope this helps,

Ced


Conversion from continuous to sampled signal? I am bit confused regarding sampling as i read different types of statements in different texts

Forexample i have a continous time signal $x=sin(t)$ defined for $t=0:10$
and i want to sample it with sampling frequency 5000hz
How can i obtain sampled signal $x[n]$ in this case

Please guide me which of following two ways is correct to define sampled signal?

$x[n]$=$sin(ts)$               where ts=1/5000

               or

$x[n]$=$sin(ts/fs)$               where ts=1/5000 and fs=5000

You typically write a sine wave as
$$x(t) = sin(2 \cdot \pi \cdot f \cdot t)$$
where t is time and f the frequency of the sine wave.

Please note that both time and frequency are physical quantities, i.e. they have values and units. The units need to cancel each other since you can't apply the sine function to a unit

If you want to sample this at $f_s = 5000 Hz$ you simply calculate the value every a time interval $T_s = 1/f_s$. You get

$$x_d[n] = sin(2 \cdot \pi \cdot f \cdot n \cdot T_s) $$

or

$$x_d[n] = sin(2 \cdot \pi \cdot n \cdot f/f_s) $$

Whether this is "correct" or not depends the ratio of frequency to sample rate and whether this meets the criteria of the sampling theorem or not.


Wavelet is a small wave, what is problem with using sin(x) of one cycle? As far as I know a $\sin(x)$ which is one cycle long can be made into a wavelet and then we can use its dilated and translated versions to represent another function, effectively a wavelet transform.

We have wavelets like Haar and others, why not have a single cycle sinusoid as a wavelet?

A 1-cycle sinewave wavelet corresponds to imaginary component basis vector for bin 1 of a rectangularly windowed DFT result, and thus has the same problems of the lowest bins near DC of a rectangular windowed FFT.

The frequency response of a rectangular windowed FFT is a Sinc function with an expanse of large ripples in the response. The main lobe of the frequency response would be very wide relative to the frequency of that wavelet's basis. And the complex conjugate image, or negative frequency response, is nearby to the lowest bin(s); and thus their side lobes can potentially cause interference with the positive frequency response, depending on phase.

e.g. the wavelet would be highly subject to interference  from noise and any other spectral content that wasn’t purely sinusoidal and exactly integer periodic relative to the wavelets length.


How to filter a white noise or random noise in Matlab Simulink? I am working on a sun-synchronous satellite Simulink model wich determines its attitude from a Sun sensor and a star sensor whom measure the sun and star direction in the satellite body reference frame. I have to introduce measure errors on both sensors using a "white noise" Simulink block, but I prefer to generate the noise through a Matlab function block using the function "rand" (or shoul'd I use the "randn" function?).

My issue is that I have to filter the noises on the measurements before sending them to a block who compute the attitude quaternion error to feed it to a PID controller blockset which commands a Reaction Wheels Assembly. Without noise the satellite flies smoothly.

I have tried some solution such as lowpass filter, Butterworth filter, mean on peaks and downs, but nothing seems to work properly.

I cannot figure out, for example, which should be the cut off frequency to set up on a lowpass filter ( I cannot figure out if the lowpass filter is the correct choice for such a problem, too).

It may be useful to have a deeper understanding of you satellite's dynamics.
What I mean is that by observing your attitude quaternion error variation DSP ( or using the motion mathematical model), you may be able to choose a convenient filtering type/tuning to get rid of most possible white noise.


How to use low pass filter in Matlab for discrete data? As see in the picture, there are high frequency noise in the signal. How to design a low-pass filter to eliminate the noise? How to apply this filter in Matlab?

Check conv function. Use conv(f,x) where f is the filter and x is the signal.

You can also implement the filter process in frequency domain:

X=fft(x);
F=fft(f);
X_LP=real(ifft(X.*F));

You can use filter too. In y = filter(b,a,x), the filter is described by numerator coefficient vector b and denominator coefficient vector a,which is different in conv where the input is the time series of the filter. But conv(f,x) and filter(f,1,x) gives the same result in your case.

EDIT

For dfilt.dffir type filter, use the code below (example page):

x = randn(100,1);    % Original signal.
b = fir1(50,.4);     % 50th-order linear-phase FIR filter.
hd = dfilt.dffir(b);    % Direct-form FIR implementation.
y = filter(hd,x);

filter can handle FIR and IIR systems, while conv only works for FIR filter.

Besides,  filter can also return the filter states, so that it can be used in subsequent calls without incurring filter transients.


Estimation of the amplitude of a frequency component in a signal I'm working on a device that produces a signal with sampling rate of 100 kHz. Every 40 ms, I receive 16000 points, from which I have to extract a frequency close to 3.55 kHz. I have developed a method that extracts the frequency of the signal very close to the Cramér-Rao lower bound, but I need to extract the amplitude of the signal at that frequency very precisely.

The first idea that occured to me was using a digital implementation of a lock-in amplifier, but the problem is that a lock-in amplifier needs a low-pass filter. And I'm trying to avoid using filters as much as I could, because filters introduce a systematic shift in my readings, while I need extremely high precission and accuracy.

What's the best amplitude estimation algorithm I should use?

Generally, the frequency estimation problem decouples from the amplitude + phase estimation problem.

As I said in the comments, you could just do:

$$
\hat{A} = \left|\sum_{n=0}^{N-1} x_n e^{-j 2\pi\hat{f} n} \right|
$$
where $\hat{f}$ is your frequency estimate.

Regarding your issues:

because the number of points isn't high enough to hit the right resolution

I'm not sure what you mean by "resolution" here? Resolving the amplitude?

second is that I need a really fast way

The suggested technique will only take $N$ complex multiplications and $N-1$ complex additions, plus whatever you need to calculate the absolute value (two complex multiplications and one addition).

That's pretty fast, and linear in $N$.  And certainly faster than an FFT (which would be $N\log_2 N$).


Time scaling of discrete-time sequences and the DTFT In the second edition of Signals and Systems by Alan Oppenheim, he discusses the DTFT of a "time-expanded" sequence that is effectively a slowed down version of the original sequence and can be thought of as the original sequence $x[n]$ with $k-1$ zeros placed between each sample. However, the DTFT for the sequence $x[an]$ where $a$ is an integer is not discussed. Why is that? Is it because information is lost when we talk about $x[an]$ for example (since it only contains samples evaluated at multiples of $a$)?
One other thing that I noticed is that unlike the CTFT, Oppenheim only discusses differentiation in the frequency domain, unlike the CTFT where both integration and differentition in the frequency domain were presented. Why is that?

The DTFT of some sequence $x(an)$ is discussed in the book you mention, but rather indirectly.

It's easy to see that $x(an)$ is a decimated sampled version of $x(n)$. This means that we take $x(n)$, we sample it with sampling period $a$, and then we decimate it (i.e. we take every $a$-th sample).

You can derive the DTFT of $x(an)$ knowing the DTFT of $x(n)$ with this in mind. If you have a copy of Oppenheim's Signals and Systems, check the section Sampling of discrete-time signals. Everything is explained there to answer your first question.


Multiple peaks in a same signal? I am working on peak detection in different signals, the signal plot looks like this:

After applying peak detection algorithm and tuning it for each signal, final output looks like this:

As you can see in figure 3 and figure 4, there are more than 2 peaks in a single point, for example in figure 3 there are two peaks from time 0.4 to 0.5 and 0.3 to 0.35 also two peaks. I am quite new in DSP, I wanted to confirm, is it ok to showcase this diagram, the algorithm is detecting the peaks correctly or should I show only one peak in figure 3 from 0.3 to 0.35 and 0.4 to 0.5 like below diagram?

Thanks!

is it ok to showcase this diagram ?

Sorry, we can't tell you this. This is really a function of your specific application, the physics behind your signal, and how you exactly define "peak".
DSP can offer you methods to suppress "near by" peaks or to de-noise multiple peaks but whether that's the right thing to do or not depends largely on what do want to do with the results and what else you already know about the signal.


Function implementation in Matlab I want to implement this function $[n,m] = a^{n+m}u[n]$ using 41×41 array, where $u[n] =1$ since all the value will be positive, and $n,m$ from 0 to 40. But it gives a wrong value.

N=40
M=40
a=0.9
x =zeros(41,41)
for col = 1:41
  for row = 1:41
    for k=0:N
      for i=0:M
        x(row,col)=a.^(k+i)*1;
      end
    end
  end
end

I believe this will do:

x=a.^((0:40)+(0:40)’)


How to analyze image quality? I have a few sequences of images, taken at fixed camera position, that address a moving object, which means the object moves towards the camera. But all the sequences of images are the cropped images of the object and all images are in the same size (grayscale images). From naked eyes, we can see the images that are closer to the camera are  better than the ones that taken at far distance. How can I show scientifically or show a proof that the images have varied quality? I also want to compare the quality between  sequences.

You are apparently in the context of no-reference, reference-free or blind image quality assessment. The topic is quite active, and I am not sure people have already a completely accepted framework for that. Multiple distortions may affect images: random noise, compression artifact, static blur, motion blur, etc. They require different metrics (benchmark data here). Some references to start with, and help you for more focused questions:

No-Reference Image Quality Assessment in the Spatial Domain, 2012
Making a “Completely Blind” Image Quality Analyzer, 2013
No-reference image quality assessment based on spatial and spectral entropies, 2014
MDID: A multiply distorted image database for image quality assessment, 2017

The first three papers come from the same group, an hint that the paint is still wet.


Cascading FM synthesis operations I'm trying to understand how to program cascading FM synthesis operations.  I can find a good deal of information out there on how to do it on various pieces of hardware, but unfortunately, there isn't enough information to be able to replicate these things in software.

I know that some details may be hardware specific, and that "replicating results" may be complex as there could be specific numbers of bits being used for calculations, or perhaps hardware that alters the results (tubes), but I'm just trying to figure out what is generally considered to be the correct way to do cascading FM synthesis operations.

Let's say that I have the concept of a free spinning oscillator object in my program which stores phase, and has a function to advance phase based on the frequency passed in:

struct SFreeSpinningOscillator {
    SFreeSpinningOscillator () : m_phase(0.0f) { }

    // advances the oscillator and returns the current value of the sine wave
    float AdvanceOscillator (float frequency, float sampleRate)
    {
        float ret = std::sinf(m_phase * c_pi * 2.0f);
        m_phase += frequency / sampleRate;
        return ret;
    }

    // note that phase isn't a true angle, but is just a percent 0..1
    float m_phase;
};

With the above I know that i can get vibrato easily enough with something like this:

SFreeSpinningOscillator a, b;

....

float modulator = a.AdvanceOscillator(10.0f, sampleRate) * 5.0f;
float carrier = b.AdvanceOscillator(noteFrequency + modulator, sampleRate);
return carrier;

The above will make it so the carrier value is a sine wave at frequency "noteFrequency" which varies it's frequency +/- 5hz on a sine wave, 10 times a second.

When talking about FM synthesis, it's more likely to be something like the below where the modulator's frequency and amplitude is based on the note's frequency:

float modulator = a.AdvanceOscillator(noteFrequency * 0.5, sampleRate) * noteFrequency;
float carrier = b.AdvanceOscillator(noteFrequency + modulator, sampleRate);
return carrier;

This is where I start to get confused though.  What if I have two or more modulators?  Are their frequencies and amplitudes ALSO based on the basic note frequency?  Or are they based on the frequency of the last modulator?

If it's based on the basic note frequency, how do the amplitudes of the modulator waves fit in?  If it's based on the modulator frequency, why even calculate the wave forms?

In other words, is the below correct? or is there some other way that it is usually or commonly done?

SFreeSpinningOscillator a, b, c, d;

....

float modulator3 = a.AdvanceOscillator(noteFrequency * 0.5, sampleRate) * noteFrequency;
float modulator2 = b.AdvanceOscillator(noteFrequency * modulator3 * 0.2, sampleRate) * noteFrequency;
float modulator = c.AdvanceOscillator(noteFrequency * modulator2 * 2.8, sampleRate) * noteFrequency;
float carrier = d.AdvanceOscillator(noteFrequency + modulator, sampleRate);
return carrier;

Here are some of the sources I've been reading and have found most helpful:

http://the-all.org/tx81z/fm_overview.html

http://www.soundonsound.com/sos/apr00/articles/synthsecrets.htm

https://www.soundonsound.com/sos/may00/articles/synth.htm

the code on the bottom looks like a triple-nested FM oscillator if i understand how the AdvanceOscillator() method works.  likely the SFreeSpinningOscillator class has a method of setting the initial phase of the oscillator.  you may get wildly different tone timbres with exactly the same code and different phases for the 4 oscillators.

if you can get your hands on a DX7 alg chart

you can see that most of their algs are sum of two chains.  alg 2 has a triple nested oscillator chain where your d is the DX7 osc #3.  but i can't tell you what the initial phases of their oscillators are.  there are so many parameters that affect the timbre, sometimes in unpredictable ways.  so doing software to emulate a DX7 will involve some tweeking of knobs


How to compare and measure the similarity between two between two signals? I'm working on different sensors data, the question is this: I have two signals (coming from two different sensors measuring two different phisical charactheristcs of the same "object").
I know for a fact that there is a certain amount of time delay between the two signals,  but beyond that, I just want to know how is the correct procedure to correlate the two signals in order to have a proper idea of the similiarity of the signals. Moreover: Can I modify the signal sens1 in order to predict (if in real time or in order to get in post-processing) the signal sens2?
Is any specific tool or package in python to easily do this type of comparison?

Thaks a lot in advance,

BR

The DTW approach is applicable to signals where there is an acceleration or deceleration between the observations of similar signals during the data capture and will show how similar or not the signals are besides the distortion from the "time-warping".

I am not sure that would be the best approach for this case given nothing is indicating there will be a change in the rate of time between the sensors (although to a second order the sampling clock in one will not be exactly matched to the sampling clock in the other so DTW could be useful if the sampling clock offset exceeded the coherence bandwidth of the capture, which is $1/T$ where $T$ is the length of the capture for which we will use correlation, which I will describe next).

This looks like an ideal application for the cross correlation function, which will show the correlation between the two waveforms for every time offset between the two. This is done by first removing the mean from each waveform, and then multiplying the two resulting zero-mean waveforms together element by element and summing the result, repeating for each possible sample shift between the waveforms. These results can be scaled by the product of the standard deviation of the two waveforms, which would normalize the result similar to what is done in the Pearson correlation coefficient.

Another similar approach is to use the Wiener-Hopf equations to solve for the equivalent "channel" between the two waveforms, which I have explained in more detail at this post which evaluated the signal received by two microphones:

Compensating Loudspeaker frequency response in an audio signal


Degraded image: Wiener Filter vs Wiener deconvolution I am studying methods for recovering images (that has been blurred and noised).

Can someone explain the difference between:

http://en.wikipedia.org/wiki/Wiener_deconvolution
http://en.wikipedia.org/wiki/Wiener_filter

I can understand that wiener deconvolution is used for Deblurring while wiener filter is used for Denoising.

Is that right? Are there any difference between the two?

Wiener deconvolution is an approach to solve the deconvolution problem that relies on the filter proposed by Wiener. The equation is the same in denoising and deblurring, except that the filter $G$ (to stick with Wikipedia's notations) that you should use is different.

To make things clear:

denoising consists in the case where the degradation kernel $H$ is the identity,
in deconvolution this kernel is equal to the blur kernel that destroyed your image (motion
blur, lens PSF...).

Hence, in denoising, the $H$ part that appears in the expression of the Wiener filter $G$ disappears (it is equal to 1), while in deconvolution it is equal to the Fourier transform of your blur kernel.

More generally, when using an inverse problem approach for these tasks, i.e. solving $\| A x - y \|^2 + \lambda g(x)$ (with $x$ your solution and $y$ the noisy observation, $g(x)$ some prior on $x$), you can always switch between denoising and deconvolution by choosing $A = Id$.


Inverse DTFT Problem Having trouble finding the inverse DTFT of $\ X(\ e^{j \omega}) = \frac{3 - \frac{1}{4} e^{-j\omega}}{1 - \frac{1}{4} e^{-2j\omega}}  $

Given the IDFT of $Xe^{j \omega}$ as :

$x(n) = \frac{1}{2\pi} \int_{-\pi}^{\pi} X(\ e^{j \omega}) \cdot e^{j \omega n} d\omega$

Tried long division as an option and reduction to partial fractions but not to any good so far.

Not sure on how to proceed.

You can use partial fractions first to expand the $X(e^{j\omega})$ as

$$X(e^{j\omega})=\frac{\tfrac{7}{4}}{1+\frac{1}{2}e^{-j\omega}}+\frac{\tfrac{5}{4}}{1-\frac{1}{2}e^{-j\omega}} $$.

Then taking the inverse DTFT you get $\tfrac{7}{4}(-\frac{1}{2})^{n}u(n)+\tfrac{5}{4}(\frac{1}{2})^{n}u(n)$


Why are two cameras mounted in paralle in stereo application? For stereo cameras on the market, two cameras are always mounted side by side and with a displacement that is perpendicular to the cameras’ optical axes. I take this setup for granted. One idea came to my mind whether this is necessary? If two cameras are not parallel and have different focal length, camera calibration can correct the difference. Why are two cameras mounted in parallel? My guess is that two cameras can have a large overlapping region. Am I correct?

When you say cameras are parallel  it only means that the center pixels are parallel.

Every pixel has a different angle, and that refers to the overlapping regions you were talking about. Actually there is at most one overlapping voxel (3d pixel kind of) for every 2 pixels!

There is more about it in epipolar geomtry , it's interesting and pretty simple!


help in understanding cosine filter I was referring to this link https://cdn.selinc.com/assets/Literature/Publications/Technical%20Papers/6059_HowMicroprocessor_Web.pdf?v=20180606-230156 . I am not very clear about the derivation of the cosine filter and how and where it comes from?

As described in the paper and observed by @Ben in the comments below, by having two filters each respectively with the sine and cosine of this filter specifically one cycle of the passband frequency of interest, we are converting what would otherwise be a simple moving average (where the coefficients are all one) to a complex bandpass filter centered at the line frequency. Complex meaning the filter has two outputs given as I and Q or 'real' and 'imaginary', and for such complex numbers the magnitude is given as $\sqrt{I^2+Q^2}$. We see this from Euler's equation relating the magnitude and phase of a complex phasor to the real and imaginary components as sine and cosine terms:
$$Ke^{j\theta} = K\cos(\theta)+jK\sin(\theta) = I+jQ$$
Which rotates the moving average filter response to the right on the frequency axis as plotted below:

The changing a moving average filter which is a low-pass with a Sinc function frequency response (Dirichlet Kernel specifically, approximates a Sinc as the number of coefficients in the filter is large) to be centered at the line frequency is an example of the process of converting a low-pass filter to a band-pass by heterodyning the coefficients. For more details on that see Different way to separate a particular frequency from a signal
The OP also asks in the comments to this answer for further insight into the filter operation with regard to the pole-zero plot. We can quickly create the pole-zero plot by inspection by knowing the result for the simple moving average filter and rotating it by angle $\omega$ consistent with the change to complex coefficients given by $e^{j\omega}$. We may already see this in the frequency response above where the aliased Sinc response of a Dirichlet Kernel has been circularly shifted to the right. To see this further, consider that a moving average filter, as an FIR filter, has only non-trivial zeros; all the poles will be at the origin for causality. For an N-tap moving average (16 in this example), the transfer function is given as:
$$H(z) = \sum_{n=0}^{N-1}z^{n}$$
The zeros are the roots of this equation, which for this case are the Nth roots of unity, omitting the zero at $z=1$. This may be clearer by knowing the roots of the polynomial $1-z^{-N}$ describing a simple differencing filter are all the N roots of unity; zeros evenly spaced around the unit circle including z=1, and we see from this post CIC Cascaded Integrator-Comb spectrum that the moving average filter can be formed by cascading this simple difference filter with an accumulator which has a pole at z=1, cancelling the zero that was at that location.

By implementing a complex filter with sine and cosine coefficients consistent with $e^{j\omega}$ where in this case $\omega = 2\pi/16 = \pi/8$, we would just rotate the pole zero diagram above by that angle to get the following result:

Essentially the filter rejects the negative frequency of the real line signal, converting the real signal that has a magnitude oscillating sinusoidally to a complex signal whose magnitude is constant for every sample, thus making it feasible to estimate the magnitude in noise free conditions for any single complex sample of the filter output (once the filter has converged from an initial transient, meaning any one sample after 16 samples for the 16-tap filter. We can average this result over more samples to improve our estimate of the magnitude in the presence of noise.
As we see in both the spectrum and pole-zero map, the filter when centered on the line frequency has excellent rejection at all the harmonics of the line frequency given the zeros are placed at all the integer harmonics. The higher sidelobes suggest sensitivity in the filter to noise at all these locations and rejection of this could be further improved by increasing it's length to include more integer cycles and windowing the coefficients with a window such as the Kaiser function. (for further details on that see:  FIR Filter Design: Window vs Parks McClellan and Least Squares)
The significance of this is demonstrated below when we observe on a dB scale the original frequency response out to half the sampling rate of 960 Hz (using the coefficients as given in the paper) with a Kaiser windowed frequency response below (which has 4x the number of coefficients at the same sampling rate). The original is in blue and has excellent rejection of harmonics but very poor rejection everywhere else, while the orange plot is the frequency response of the cosine and sine filters with the additional Kaiser windowing.  This filter would maintain a much better estimate of the magnitude in the presence of non-harmonic noise, as well as provide excellent rejection of harmonics.


Converting Time Shift to Phase Shift I have just started taking signal & systems lessons and here is my question:

If we say that x0(t) = A.cos(w0t). (a cosine signal with zero phase shift, w0 radian frequency) And if we would like to time-shift this signal by for example t1.

The resulting phase shift should be found by -w0*t1 = phi. (phi denoting phase shift).

By this logic, if I wanted to find the phase-shift for this signal: x(t) = sin(pi/4*t-pi/4)

What would the result be? Can I apply the first formula for phase for already phase shifted signals? Or am I missing something very big here? :)

Thank you very much for any kind of help.

Indeed, you can apply the same technique:

As you said, given a sine signal with time shift $t_1$, which is $x(t)=\sin(w(t-t_1))$, this becomes $x(t)=\sin(w t - wt_1)=\sin(wt+\phi)$ and $\phi=-wt_1$.

Now, if you encounter a signal $x(t)=\sin(\frac{\pi}{4} t - \frac{\pi}{4})$ you identify $w=\frac{\pi}{4}$ and $\phi=-\frac{\pi}{4}$. What remains is to solve $\phi=-w t_1$ for $t_1$ and you have it.


sinus generation and it's fft in radians I'm trying to generate sine signal and display it's fft in radians, byt I meet several issues that I can't understand. Please find my script below:

clear all; close all; clc;
fs = 10*pi;
ts = 1/fs;
t = 0:ts:2*pi;
sine = cos(10*t);
plot(t/pi,sine)
sinefft = fft(sine);
subplot(2,1,1)
plot(abs(sinefft))
subplot(2,1,2)
plot(abs(fftshift(sinefft)))

Results

1. Why on first plot fft shows that cosine has a frequency 11 radians?

2. How to shift fft to achieve frequency domain around 0?

I will be thankful for your explanations.

Maks

Why on first plot fft shows that cosine has a frequency 11 radians?

The array index number of the cosine peak is 11.  The frequency is different.

Matlab indexes from 1, so the FFT bin number will be 10 (as we have a 0 bin in FFTs).

The frequency of that bin will be
$$
\frac{2 \pi (11 - 1) }{N } = \frac{2 \pi 10}{ 198} = \frac{10\pi}{99}
$$
assuming your FFT length was 198.

Note that this will be a normalized frequency between 0 and $2\pi$ (or, if you prefer, $-\pi$ and +$\pi$).

  How to shift fft to achieve frequency domain around 0?

Your use of fftshift does this in the second plot.  However, you need to set the $x$ axis by doing something like:

N = length(sinefft);
delta_omega = 2*pi/N; % rad/sec between FFT bins
omega = -pi:delta_omega:(pi - delta_omega)
plot(omega, abs(fftshift(sinefft)));

Minor Correction I previously had this frequency variable as f, it's actually omega.

  How to calculate now what frequency in Hz has this signal?

So the normalized frequency in radian per second is omega.  That ranges from $-\pi$ to $+\pi$ (modulo some end effects).

The normalized frequency in "Hz" is f. This will be omega/2/pi. This ranges from -0.5 to +0.5.

The frequency in terms of the sampling frequency will be f_true = omega / pi / 2 * fs (where fs is the sampling frequency).


Viewpoint dependant deformation of objects in photos I have a question concerning the deformation of objects in photos. Which deformations can an object depicted on a photo undergo through viewpoint and condition changes?

Some deformations I have in mind are:

rotation
scale
shearing
translation
illumination
occlusion

Which deformations did I forget?

Best!

Since you have listed "illumination", I am assuming that "condition changes" include changes in lighting.  So I would add cast shadows, which will change if you move the light source.

Also, all the geometric transformations you have listed are affine. But if you move the camera perspective changes as well. So you may want to consider the projective transformation.


Detection markers in images I have a set of images with red, round labels on it. What I want is, to get pixel location of each of the label. One way, of course, would be to search for red color pixels. But the problem is, sometimes light reflexes from the marker and changes it's color. But you could still detect round "redish" circle. It there any efficient way to get this?

You can always use normalised cross correlation with a structuring element like a disk or ellipsis to account for different viewing angles and marker distances. This is basically template matching, where a target that "looks like" the template is searched for in a larger image.

If more computational power is available, you could go down the route of detecting circles via the Hough Transform.

In this case, the Hough transform is used to identify some strong features that appear as points in some space and then, based on assumptions about their locations, the position of the circles is inferred. Similar to those techniques, are others that identify generic features off of an image (e.g. corners) and then try to cross-reference the position of the identified features with the relative positions identified off of a template.

A compromise might be to use similar methods as used to detect shotCodes. Now, I did try to locate that paper by the original authors of shotCodes where the detection process was outlined in detail but I cannot seem to find it. Instead, here is one that talks about detection (and the (more recent) associated thesis). This was quick enough to be implemented on early smartphones. Perhaps it offers a compromise in terms of ease of detection and position accuracy off of images.

Hope this helps.


Given an FIR filter coefficients, how can I find coefficients of its equalizer filter? I have an N-tap direct form complex FIR filter with known coefficients [c0 c1 .. cn]. I want to find the coefficients of its equalizer filter which when cascaded with the original filter gives constant magnitude response for all frequencies. How can I do that? Thanks.

Assuming that you also want to equalize the filter's phase response (not only its magnitude response), you need an equalizer with a transfer function $E(z)$ that is the inverse of the FIR filter's transfer function $H(z)$:

$$E(z)=\frac{1}{H(z)}=\frac{1}{c_0+c_1z^{-1}+\ldots+c_{N-1}z^{-(N-1)}}\tag{1}$$

This is an all-pole filter which can be implemented by the following recursion:

$$y[n]=\frac{1}{c_0}\left(x[n]-c_1y[n-1]-\ldots-c_{N-1}y[n-N+1]\right)\tag{2}$$

where $x[n]$ is the input sequence, and $y[n]$ is the output sequence.

The problem with implementing Eq. $(2)$ is that the corresponding filter can be unstable. This depends on the locations of the zeros of the FIR filter's transfer function $H(z)$. If $H(z)$ has zeros on or outside the unit circle of the complex plane (i.e., zeros with magnitudes equal to or greater than $1$), then $E(z)$ has poles at those locations, and poles on or outside the unit circle cause instability of the equalizer. This means that if $H(z)$ has zeros with magnitude equal to or greater than $1$ (i.e. if $H(z)$ is not strictly minimum phase), then there exists no causal and stable equalizer.

In practice you can design an approximate equalization filter by choosing the filter coefficients of the equalizer such that the error

$$\epsilon=\sum_i\left|H(e^{j\omega_i})E(e^{j\omega_i})-e^{-j\omega_i D}\right|^2$$

is minimized, where $D$ is some appropriately chosen delay, and the sum is computed on a dense frequency grid. If the equalizer $E(e^{j\omega})$ is chosen to be an FIR filter, then the resulting optimization problem can be solved by solving a set of linear equations.


How to estimate for an exponential in WGN Question:

Formulating a Maximum likelihood estimator:

So, the likelihood will be

$p(y;\mathbf{h}) = \frac{1}{{(2 \pi \sigma_\eta^2)}^{T/2}} \exp{(-(y - y_0(t))^2)/ 2\sigma_\eta^2}$.

Then, I need to differentiate w.r.t the unknowns and equate to zero.

This is a nonlinear equation in $x$ and cannot be solved directly. Newton-Raphson is a method but it works only for very close initial guess to $x$ and $d$. So, Was thinking how to apply Expectation - Maximization.

I don't know if my approach is correct or not.

Maximum $d$ you may get is a total number $N$ of $y_n$'s you have - 1. Larger $d$'s yield curvatures that can't be represented by N points.

Solve a matrix equation $A\vec{t}=\vec{y}$ where $\vec{y}$ is a measurements vector (size $N \times 1$) and $A$ is a normalized version of matrix $X$ (size $N \times N$), where each column $A_j$ is calculated as:

$A_{ij} = \frac {X_{ij}} {\sum \limits_j (X_{ij})^2}$

$X = \left [ \begin{matrix}
x_0^0 & \cdots &x_0^{N-1} \\
\vdots & \ddots & \vdots \\
x_{N-1}^0 & \cdots & x_{N-1}^{N-1}
\end{matrix} \right ]$

d you are interested in is a position of largest coefficient in a vector $\vec{t}$ (size $N \times 1$ starting from zero). Your approximation would be $y = t_d \cdot x^d$.

You can achieve better results if your solver finds solution optimal in $L^1$ sense.

This is a sort of exhaustive search, yes.


Peak to Average Power Ratio of an exponentially increasing signal I am having an exponentially increasing sinusoidal signal.
After sampling, I am dividing the signal into frames of length 128 and doing windowing operation with a hanning window of length 128. Then the FFT of this windowed sequence is found out frame by frame.
The average power and peak power of the windowed signal in each frame is found out and Peak to average power ratio (PAPR) is calculated.
Will the PAPR increase with each frame?

This feels a lot like there's quite a bit of homework you should do about this, but let's give you a start:

Write down your signal explicitly as

$$s(t) = p(t) \cdot e^{\lambda t}\text,$$

with $p$ being your periodic oscillation. You sample that, and divide it into vectors. Each vector has $N=128$ samples. Let's write it this:

\begin{align}
s_k[n] &= s\left(\frac{n+N\cdot k}{f_\text{sample}}\right),& n = 0,\ldots, N-1\\
&= p\left(\frac{n+N\cdot k}{f_\text{sample}}\right) \cdot e^{\lambda\frac{n+N\cdot k}{f_\text{sample}}}
\end{align}

where $k$ gives us the frame number.

Your question is how things develop when you increase the frame number!

So, consider the quotient of the next frame and the current frame:

\begin{align}
\frac{s_{k+1}[n]}{s_k[n]}
&= \frac{p\left(\frac{n+N\cdot (k+1)}{f_\text{sample}}\right)}{p\left(\frac{n+N\cdot k}{f_\text{sample}}\right)} \cdot \frac{e^{\lambda\frac{n+N\cdot (k+1)}{f_\text{sample}}}}{e^{\lambda\frac{n+N\cdot k}{f_\text{sample}}}}\\
&= \frac{p\left(\frac{n+N\cdot (k+1)}{f_\text{sample}}\right)}{p\left(\frac{n+N\cdot k}{f_\text{sample}}\right)} \cdot e^{\lambda\frac{n+N\cdot (k+1)}{f_\text{sample}}-\lambda\frac{n+N\cdot k}{f_\text{sample}}}\\
&= \frac{p\left(\frac{n+N\cdot (k+1)}{f_\text{sample}}\right)}{p\left(\frac{n+N\cdot k}{f_\text{sample}}\right)} \cdot e^{\lambda\frac{N}{f_\text{sample}}}\\
\end{align}

Now, it'd be a nice simplification if we can just restrict the oscillations' period to something that is an divisor of $N$. That especially means that $p()$ becomes $N$-periodic (under the given sampling rate).

\begin{align}
&= 1 \cdot e^{\lambda\frac{N}{f_\text{sample}}}\\
&=\text{const.}
\end{align}

In other words, the next frame is identical to the previous frame, aside from a constant factor.

Since scaling with a constant factor doesn't change the PAPR at all, this was easy to answer for oscillations that have a frequency that divides the sampling rate, only. The rest isn't really unusual (it's leakage, essentially), but the math is slightly more annoying, so I'll leave it to you.


How the noise in time series can be reduced? I have asked the question in the following link: Are time series data always contain noise?
I wonder if we could use variational methods to denoise the time series data just as in computer vision?
Thanks in advance...

Yes indeed, there is a huge body of works on the topic of variational methods for signal restoration (more general than denoising), one example being total variation denoising, or proximal methods.

If you can be more specific, you could get more accurate and precise answers.


Instantaneous magnitude I was reading this thread: How to get Instantaneous Magnitude for a Instantaneous Frequency From FFT?

I have basically the same question, but I'm curious about more detail about the Hann window kernel. ederwander said that once you know the frequency offset from the nearest bin then the instantaneous magnitude is a the Hann kernel of your offset times FFT length, times magnitude of the bin you're offset for.

Where the Hann window kernel $W(k) = 0.5 (\operatorname{sinc}(k(M/N)) / (1 - (k(M/N))^2)$

I'm just not sure why that kernel evaluated at offset*N gives the proper scaling factor for adjusting the magnitude. I couldn't even find any other information on the Hann kernel or how it is derived.

Thanks!!

The Hann Window Kernel is the FFT of the Hann Window.

Offset is defined strangely, and unnecessarily. calculating InstFreqs and then dividing it by Sampling Frequency Fs is just throwing away the work of attaining it.

Bin Resolution = Fs / N
//InstFreqs will start at 0Hz and peak at Fs Hz (or half that??)
InstFreqs(i) = i * Fs / N

So in the referenced answer:

Offset = abs(InstFreqs(i) / (Fs/N) - i)

but if I'm correct that: InstFreqs(i) = i * Fs / N then:

Offset = abs(InstFreqs(i) / (Fs/N) - i) = abs( i - i ) = 0

So that line didn't make sense. Too bad there was no comment hinting at what it was supposed to mean.

index is put into Wk it must span the length of the Window, and not the length of the FFT. So let's start by changing that. N * Offset will (likely) exceed the length of M (window size)

index = floor(abs(Offset) * N) + 1
InstMag = Mag(i) * Wk(index)

Unless of course I misinterpreted what InstFreqs(i) is, which it seems I have.


use prony method for damped exponentials let us consider represented of damping exponential model by prony method,there is source code

y=zeros(1,N);
for i=1:N
    y(i)=x(800*i);
end
d=zeros(1,N/2);
for i=1:N/2
    d(i)=y(i+N/2);
end
D=zeros(N/2,N/2);
for i=1:N/2
    for j=N/2:-1:1
        D(i,-j+N/2+1)=y(i+j-1);
    end
end
a=pinv(D)*d';

muhat=roots([1,-a']);
U=zeros(N,N/2);
for i=1:N
    for j=1:N/2
        U(i,j)=muhat(j,1)^(i-1);
    end
end
C=pinv(U)*y';

and equation of model is following with there solving  procedures :

with their  description :

and solving strategies :

Find Roots of charactreristic polynomial formed from the linear prediction coefficients

Solve the original set of linear equations to yield the estimates of the exponential amplitude and sinusoidal phase

i want following thing: we have given signal and i want input  be it's length and  L,or function should have following form

function [amplitudes,damping_factor,phase,frequency]=prony(y,n,L)

%n-length(y)
%L-number of complex exponential,

how can i continue?how to  change given code in my case?thanks in advance

you have C and U, and C can be written as C_i = (A_i/2)*exp(j*phase)
U_i can be written as exp(damp_i + j*2*pi*f_i)*T
its should be easy to extract the info from that. T is sampling period.

%to get damping factor i and frequency i use something like this
%assume damping factor at sample 1 is 0.3, frequency at 1 is 20
%and out sampling frequency is 500 Hz
damp1 = [0.3 0.4 0.5];
f1 = [20 30 40];
T = 1/500;
%assume this is your ui vector from th ealgorithm
ui = exp((damp1 + j*2*pi*f1)*T);
%get the damping and frequencies
lv1 = log(ui);
r_damp1 = real(lv1)/T;
r_f1 =imag(lv1)/(T*2*pi);

%assume this is your Ci vector from algorithm
Ai = [0.1 0.2 0.3];
phasei=[pi/8 pi/4 pi/2];
Ci = (Ai/2).*exp(j*phasei);
%extract phase and amplitude
r_Amp = abs(Ci)*2;
r_phase = angle(Ci);


Length of DFT and Number of Filters This plot represents the magnitude of the frequency response of certain DFT sub-band filters

I have some difficulties in understanding correctly the DFT. How can i know, only by observing the plot , the number of filters and the lenght of the DFT?

Also several times the DFT was presented to me like a parallel bank of sub-band filters and i cannot understand the full meaning of that phrase.

Thanks by the help.

One first suggestion for you is to implement DFT filter banks in some language (Octave, Python, Matlab) and play with them to get some insight. Then,
let me provide a information source, and a few hints.  In  The DFT Filter Bank, you can find:

  In this section, we will show how the DFT can be computed exactly from
  a bank of $ N$ FIR bandpass filters, where each bandpass filter is
  implemented as a demodulator followed by a lowpass filter

The above website is really worth reading. You can see for the plot that you have two responses, one being obtained via a shift in frequency for the other. Hence you know that, from one of the Fourier modulation properties, that a shift in frequency corresponds to a modulation in the time domain. Let us now look at the low-pass filter only (in blue). A running low-pass filter has a frequency response related to the sinc function. More precisely, it is called a periodic or an aliased sinc, with frequency response:

\begin{eqnarray*}
H(e^{j\omega}) &=& \frac{1-e^{-j\omega N}}{1-e^{-j\omega }}
= \frac{e^{-j\omega N/2}}{e^{-j\omega /2}}
\frac{\sin(\omega N/2)}{\sin(\omega /2)}\\
&= &
Ne^{-j\omega(N-1)/2} \hbox{asinc}_N(\omega ).
\end{eqnarray*}
where $N$ is the length of the running sum. It can be inferred  from the number of ripples, corresponding to a number of zeroes of the aliased sinc.
The number of filters is potentially trickier, because of oversampled filter banks, but this may be outside the scope of your lecture. Look at how many shifted (modulated) versions of the low-pass can cover evenly the spectrum.

Then again, do not hesitate to verify via simulations.


How to separate Even and odd part (in-phase bit streams I and quadrature-bit streams Q) of the signal? I am working on communication system in simulink and wondering, how can i separate the Even and odd part of the signal.
In my Transmitter block, I have :
$$
\boxed{\text{Binary Data signal}}{\longrightarrow}\boxed{\text{Serial to parallel}}^{\nearrow\boxed{\text{Even Part}}{\longrightarrow}\boxed{\text{Multiply PN code}}{\longrightarrow}\text{signal}}_{\searrow\boxed{\text{Odd Part}}{\longrightarrow}\boxed{\text{Multiply PN code}}{\longrightarrow}\text{signal}}
$$
Not sure how to separate those two part from signal ? Also can i use below equations for channels I and Q.
$$
I(t) = \sqrt{\frac{2E_s}{T} \cos(\omega t)\cos(45+n)}
$$
$$
Q(t) = \sqrt{\frac{2E_s}{T} \sin(\omega t) \sin(45+n)}
$$

The OP is mis-understanding how this QPSK signaling scheme operates.
There is a stream of input bits $b_0, b_1, b_2, b_3, \ldots,$ that is coming serially at the rate of one bit every $T$ seconds, so $b_0$ during $[0,T)$, $b_1$ during $[T,2T)$, and more generally, $b_n$ during $[nT, (n+1)T)$.. So, there is a clock signal at $\frac  1T$ Hz at the input.
The first thing that happens is the serial-to-parallel conversion where the serial bit stream gets converted into a stream of dibits (two bits in parallel), one bit on each of two individual wires and each bit lasting for $2T$ seconds (that is, the dibit lasts for $2T$ seconds.  Thus, we have
$$(b_0, b_1) \to \left[\begin{matrix}b_0\\b_1\end{matrix}\right], \quad (b_2, b_3) \to \left[\begin{matrix}b_2\\b_3\end{matrix}\right], \cdots \cdots
(b_{2n}, b_{2n+1}) \to \left[\begin{matrix}b_{2n}\\b_{2n+1}\end{matrix}\right], \cdots$$
So immediately we need to understand that

The output of the serial-to-parallel converter is clocked at with a clock at $\frac{1}{2T}$ Hz, and this clock is synchronized with the input clock of the serial to parallel converter,
There is a delay of $2T$ seconds in that $\left[\begin{matrix}b_0\\b_1\end{matrix}\right]$ is the output of the serial-to-parallel converter during $[2T,4T)$,  $\left[\begin{matrix}b_2\\b_3\end{matrix}\right]$ is the output during $[4T,6T)$, and so on,
The "even bits" have been sorted into appearing on the upper wire and the "odd bits" on the lower of the two wires, just as is envisioned in the OP's diagram. These are the Inphase and Quadrature branches of the transmitter.

The dibit enters into the direct-sequence modulator which says multiplication but this is if the bits ($0$ or $1$) have been converted into real numbers $+1$ and $-1$ respectively, and ditto the PN sequence. The logical operation with $0$s and $1$ is the XOR operation: the bit $b$ lasting $2T$ seconds is replaced by a serial stream of $N$ bits $(b\oplus \text{PN}_0, b\oplus \text{PN}_1, \cdots, b\oplus \text{PN}_{N-1})$ which last for a total of $2T$ seconds.  Thus the PN sequence generators (usually different ones in the I and Q branches) are clocked at $\frac{N}{2T}$ Hz, and this clock is also synchronous with all the other clocks.
Everything is digital up to this point. To get the actual transmitted signal, the $0$s and $1$s are mapped to $\pm 1$, each rectangular pulse (of duration $2T/N$ seconds in the direct-sequence QPSK system) is replaced by the desired baseband pulse shape (e.g. root raised cosine) to get the complex baseband representation $I(t)+jQ(t)$ of the QPSK signal. The actual RF signals in the I and Q branches
of the transmitter are $I(t)\cos(2\pi f_c t)$ and $-Q(t)\sin(2\pi f_c t)$ where $f_c \gg \frac{1}{2NT}$ is the RF carrier frequency. Thus the actual RF QPSK
signal is of the form
$$s(t) = I(t)\cos(2\pi f_c t)-Q(t)\sin(2\pi f_c t) = \text{Re}\big((I(t)+jQ(t))\exp(j2\pi f_ct)\big)$$ and results from modulating a (spread) dibit onto an RF carrier. Dibits consist of two bits; all Americans know that make two bits make a quarter, which explains why they think that the Q in QPSK stands for quaternary instead of quadrature or quadriphase.....


My impulse response does not tend to zero I am doing an IFFT of frequency response data achieved with Simulation tools.

When I plot my impulse response it looks wrong as the response does not tend to zero as it should. Instead there seems to be a rise in pressure towards the end of the signal.

My script looks like this:

close all;
clear all;
A = load('PvaluesAtReceiver.txt');
f = A(:,1);
Preal = A(:,2);
Pimag = A(:,3);
Prealpad = cat(1,Preal,zeros(48000-180,1));
Pimagpad = cat(1,Pimag,zeros(48000-180,1));
fpad = cat(1,f,zeros(48000-180,1));
enter preformatted text here
j = sqrt(-1);

Prealf = Prealpad(88:177);
Pimagf = Pimagpad(88:177);
ff = fpad(88:177);
Prealfpad = cat(1,Prealf,zeros(48000-90,1));
Pimagfpad = cat(1,Pimagf,zeros(48000-90,1));
ffpad = cat(1,ff,zeros(48000-90,1));
X_ff = Prealfpad + j*Pimagfpad;
X_ff = [zeros(ff(1),1); X_ff];
IRf = ifftshift(X_ff);
x_tf = ifft(IRf);
figure,plot(real(x_tf)), title('Impulse response filtered');

This produced the following plot:

Does anyone know why this doesn't tend to zero??

This is the frequency domain i want to transform.

Note that the DFT transforms a periodic sequence into a periodic sequence. Hence, the result of your IFFT is intrinsically periodic and you should view it as such. The parts you are seeing at the end of it are thus just as well behind the start as they are before it.

Hence, this part may as well be some form of pre-ringing, stemming from a nonzero group delay of your filter.

Note also that the FFT implicitly periodifies your spectrum. It looks a little cut, so you may be missing some essential parts.


DFT detemine frequencies I am just begining learning about DFT and I am bit unsure on what is happening mathematically during DFT.

I've sampled a signal , a sine wave sin(1000*2*pi*t). And performed DFT to calculate the Ck-coefficients, by plotting the Ck coefficients with different k values.

The frequency resolution is 8000/205
where 205 is the amount of datapoint/samples and 8000 sample frequency.

I can by using the plot read, where |Ck| is highest and then using my frequency resolution, detemine the frequency.

I see that the amount  of samples have alot to say about determining the frequency.
If i use a lower amount of samples, or does not use all the samples, will the estimated frequency either be lower, or higher. Why does is do that?

ex.

k := Round[1/((8000/100)/1000)]
freqres := 8000/100 // N
estimated freq = 468.293

k := Round[1/((8000/205)/1000)]
freqres := 8000/205// N
estimated freq = 992

k := Round[1/((8000/500)/1000)]
freqres := 8000/500
estimated freq = 992

The FFT get a discrete spectrum, as a result do you have a finite number of frequency or bins, this means that FFT resolves a time waveform into its sinusoidal components, the frequency resolution is explicitly linked to the FFT size and the sample rate, the resolution frequency is:

$$ \frac{Fs}{ N}$$

Where Fs is your sample rate and N is the size of your FFT, The larger the size of the FFT more spectral component you have and it makes its results more accurate, this tells me that to get lower frequency resolution you need a large size FFT (2048, 4096).

Follow your example 8000/205 this tells me that for each FFT bin returned you may be losing precision in order 39.02439024390244Hz


Are there any recursive online max/min filters for time-series Are there any online recursive filters that can approximate local, time-varying minimum and maximum values of a time series?

I think you are looking for a "peak" detectors. They are not filters in the conventional sense. Typically the output is increased when the incoming sample is either above a certain threshold or above the current output and decreased if it's smaller. The increase and decrease are often just simple exponential with time constants determined by your application. The increase time constant is typically called "attack" time and the decrease time constant is "decay" time.


How to find the inverse Z transform of this function in z domain? i have a function
$$F(z)=\frac{z-0.4}{z^2+z+2}$$ i need to find the inverse z transform of it , i have tried it with residues but the roots are too much ugly and it involves lots of messy calculations , is there any other way to proceed ?

I think you are using the partial fraction expansion in wrong way. You have to know for $z$ transform we expand it as a function of $z^{-1}$ and not $z$. for your special problem factor out  $z^2$ from denominator then divide the nominator by that.

Then you could factor your polynomial over a real field or complex field. This means after factorization over complex field you may obtain polynomials with complex coefficient but factorization over a real field you obtain polynomials with real coefficients. so for your case the denominator is irreducible over real field but we could factor it over a complex field, which lead to a pair of conjugate poles (roots). Now you could easily obtain the inverse using the inverse of a single pole (don't worry about the complex numbers, after summing response of 2 conjugate pole the complex parts will cancel out each other and only a real part remain).
Also I have to say, for your case you will reach 2 answer considering the region of convergence (RoC) in $z$ space.


RTL-SDR units, $\textrm{dBm}$ I am using the Python rtlsdr library to get complex samples from my RTL-SDR USB stick.  I am getting power by squaring the real and imaginary components of these values, and adding them.

What unit is the result in?
If not $\textrm{dBm}$, how can I get to $\textrm{dBm}$?

It's very hard to associate the numbers you're calculating with the actual signal power at the antenna. The SDR will filter, amplify, and AGC the signal, and then it'll convert it to digital. At best, you can say that the power you're calculating is somewhat related to the true signal's power.

The power you're calculating is in linear, not logarithmic, scale. To convert to decibels, use the usual formula. However, be aware that the reference inherent when using decibels is not well specified in this case.


Separating waves of very close wavelengths I'm new to wave problems, so excuse my simplified words.
I have a signal consists of high range of frequencies . I applied a Band-pass filter to this signal and I have been able to separate it to bands of specific frequencies.

The problem that I want to solve is as follows:
Each band have a lot of individual waves of VERY, VERY close wavelengths. For instance 4.08 Hz and 4.09 Hz. even we may have 4.001 Hz and 4.002 Hz as well. The signal data is limited and I can't push the filter to any better resolution.

The question now, How I can separate the individual waves in each band? And if possible please refer me to any books which deal with such problem.

Thanks in advance,

You have an incorrect approach here. Instead of trying to design small range filters, you should just do frequency-domain analysis.

Say you have a signal $x$ consisting of $N$ samples, sampled at the rate of $f_s$ Hz. For simplicity, let's take an example with $N = 1000$ and $f_s$ = 1000 Hz. The Nyquist frequency is half of the sampling frequency, which is $500$ Hz in this case. This means that the signal you are sampling should not have frequencies above that (otherwise there would be aliasing).

Now, if you want to know what frequencies you have in your signal, compute the discrete Fourier transform (DFT). DFT is defined as

$
X(k) = \sum_{n=0}^{N-1} x(n)e^{-2\pi k\frac{n}{N}} \quad 0 < k < N-1
$

From the DFT, you get two things. First, $|X|$ is called the magnitude spectrum, which is of your interest. The second thing would be the phase spectrum, given by $arg(X)$. DFT is symmetric from the middle for real-valued inputs, so the frequency resolution you get in our example is $1$ Hz. Now, say that you would also want to know the half frequencies. Looking at the definition of DFT, you can see it is $N$-point. If you want, you can do it likewise for $M$-points. Thus, to get the half frequencies as well, take $M=2000$.

In the case of your question, what you need is a really long DFT.

Other things to note: do not compute DFT from its definition. Use a fast Fourier transform (FFT). It's faster, and more precise.


Even and Odd Filters and Quadrature I am a mathematician trying to understand some signal processing terms. My first question regards the terms 'even' and 'odd' filters. What does this mean? Does it reference the fact that upon a Fourier transform, there is a cosine and a sine function, which in turn are even and odd? If not, a link to a reliable source would be appreciated.

Second question: What is a quadrature? I have heard the term used for the phase shift between two Fourier transforms, the method for dividing a frame into four parts, and taking a gradient-like subtraction from vertical, horizontal, and both diagonal directions. Is it all of these at once?
Also, a link to a reliable source would be great.

An even function is symmetric:

$$f(x) = f(-x)$$

An odd function is anti-symmetric:

$$-f(x) = f(-x)$$

A signal in quadrature with another has all its Fourier series components phase-shifted by 90 degrees. The two signals are thus orthogonal.

To get the quadrature signal we apply the Hilbert transform (or Riesz transform in more than 1 dimension)

$$g(x) = \pm \mathcal{H}[f](x)$$

A cosine wave is even and symmetric. A sine wave is odd and anti-symmetric. Therefore if you apply the Hilbert transform to an even signal (cosine Fourier series components) you get an odd signal (sine Fourier series components). If you apply the Hilbert transform again you get the negative original signal. Apply again you get the negative of the original Hilbert transformed signal. Apply again you get the original signal. This is where the "quad" part comes in.

This is useful in signal analysis. Since the Hilbert transform only phase-shifts the Fourier series components, the energy of the signal remains constant. We can thus reconstruct the original signal (e.g. wavelet denoising). Also, given an even filter that responds to peaks and troughs, we can create an odd filter that responds to edges. The pair of filters is called a quadrature filter pair and allows the analysis of signal features. See the analytic signal for more details.


Explanation of Spectral Line Splitting when AR process is overmodeled In the book Statistical Digital Signal Processing and Modeling by Monson Hayes, it is shown (in section 8.5.1) that when an AR(2) process described by the following difference equation
$$x(n) = 0.9x(n-2) + w(n)$$

is modeled with two different values of model orders $p = 4$ and $p = 12$, a single peak appears in the power spectrum for the former case and two peaks occur for the latter one. How do I explain this?
What is the relation between the model order and the number of peaks. Please forgive my lack of proper fundamentals, but I also wish to know why I am supposed to expect only a single peak to occur and not more than one?

The equations involved in estimating model parameters $a_{p}(k)$ and $b(0)$ and power spectrum of the AR process are given below. $r_{x}(k)$ refers to the autocorrelation sequence.

\begin{equation}
    \hat{P}_{AR}(e^{j\omega}) = \frac{|\hat{b}(0)|^{2}}{\left|1 + \sum_{k=1}^{p}\hat{a}_{p}(k)e^{-jk\omega}\right|^{2}}
\end{equation}

\begin{equation}
    \begin{bmatrix}
    r_{x}(0) & r_{x}(1) & r_{x}(2) & \dots & r_{x}(p)\\
    r_{x}(1) & r_{x}(0) & r_{x}(1) & \dots & r_{x}(p-1)\\
    r_{x}(1) & r_{x}(0) & r_{x}(1) & \dots & r_{x}(p-2)\\
    \vdots & \vdots & \vdots & \vdots & \vdots\\
    r_{x}(p) & r_{x}(p-1) & r_{x}(p-2) & \dots & r_{x}(0)\\
    \end{bmatrix}
    \begin{bmatrix}
    1\\
    a_{p}(1) \\
    a_{p}(2) \\
    \vdots \\
    a_{p}(p)
    \end{bmatrix}
    =
    \epsilon_{p}
    \begin{bmatrix}
    1\\
    0 \\
    0 \\
    \vdots \\
    0
    \end{bmatrix}
\end{equation}

\begin{equation}
    |b(0)|^{2} = \epsilon_{p} = r_{x}(0) + \sum_{k=1}^{p}a_{p}(k)r_{x}^{*}(k)
\end{equation}

One way of doing power spectrum estimation of random processes actually requires finding a set of ARMA model parameters that best approximates the given observed random process $x[n]$.

In your case, you selected an all-pole signal model with $p$ poles as the roots of the denominator polynomial $A_p(z)$ with coefficients $a_p[k]$, and a $b(0)$ related with power scaling of the process. Each pole represents a peak in the spectrum for $0 \leq \omega \leq \pi$.

So, PSD estimation is reduced to finding (estimating) all-pole model parameters $a_p[k]$ and $b[0]$, from the given data by first estimating the auto-correlation sequence values $r_x[k]$ for $k=0,1,2,...,p$  , and then solving a set of linear equations (Prony or Yule-Walker etc) to find those $a_p[k]$ and $\epsilon_p$.

Note that there are a number of approaches in estimating the ACS $r_x[k]$ of $x[n]$, such as Auto-Correlation, Covariance, Burg's, Modified Covariance, etc...

Certain of these techniques (such as the autocorrelation and Burg's method) produce an artifact known as spectral splitting when they are forced to fit a high order model to an actually low order data; aka a model-data mismatch problem. There is not much you can do to prevent it. I don't know an analytic derivation of the relation between the chosen model order $p$ and the number of splitted peaks you will observe. That's also dependent on the observed noise instance being processed. So that number of peaks is also random.

However, it's stated that modified covariance method of ACS estimation can overcome the spectral splitting problem.

Furthermore, if you continue reading, you will see in the section named Selecting the Model Order (p. 445, Eq.8.23) that by trying to make a good guess of model order you may reduce the chances of generating that line splitting problem.


Plotting impulse functions on the same time axis with other functions in MATLAB I am trying to plot the following function:

I have to plot another function along with it all on the same plot. Got the other function working, but I can't seem to get this one working. The code below is my attempt at plotting the function:

t = 0:1:10;

for t=>1 && t<=1

    n = 0:1;

end

plot(t,n)

Need help

The use of logical operators is not correct. A more generic verison is:

timeStart = -3; timeStop = 2; nSample = 1000;
time = linspace(timeStart,timeStop,nSample);
timeUp = 0; timeDown = 1;
data = zeros(size(time));
data((time>=timeUp) & (time<=timeDown)) = 1;
plot(time,data,'x-')


How to Estimate the Noise of an Image / Estimation of the Noise in an Image? I want to know how to estimate the noise of an image?
Also is this represent the camera's noise?
Does anyone know how to estimate the noise of an image or camera?
Update:
Especially, I want to know about the A, B, C from as following.
How Can I plot these A, B,C curve?
Also I want to know what is C. What does C mean?

There are various methods to estimate the noise of a signal (or adapt a threshold to be used later in PCA for example).

There are adaptive methods and statistical methods. Statistical methods employ a statistical model of the noise and estimate from the data, while adaptive methods iteratively filter the data until a certain threshold of (reduced signal) accuracy is reached. These are just generalities of the most used approaches.

The following references deal (mostly) with images from cameras:

Noise Estimation from a Single Image
Automatic Estimation and Removal of Noise from a Single Image
A noise-estimation algorithm for highly non-stationary environments


Numerical higher order derivatives and time axis I have a rather elementary question. Suppose we wish to study even-derivatives of an instrumental signal say second fourth and sixth derivatives and plot it as a function of time. With each successive differentiation we lose the data points. If the initial number of data points were 1000, the second derivative will have 9998 points, the fourth derivative will have 9996 points and so on. This can be seen with the diff command in Matlab.

How should we align these derivatives with the time axis which has 1000 points for plotting purposes? The typical way is to leave the first two points of the time axis for the second derivative and first four points for the fourth derivative. Is there a mathematical justification for this way of alignment since a true derivative is defined at a point rather than this numerical approximation? In some chemistry texts, authors suggest to average the values of the x-axis and plot the derivatives with respect to that.

Thanks.

When considering finite differences there are three of them: Forward, backward and central.

Neither of them is "shorter" than the main signal. The zero is not exactly zero, it should be $\emptyset$. That is, "I don't know".

In backward difference, $y[n] = x[n] - x[n-1]$, but at $n=0$, $n-1$ is not even defined and the same applies for $y[0]$. We usually do $y[0]=0$, (which implies that $x[-1]=x[0]$).

Notice here that we wrote $y[n]$. In finite differences, in the discrete world, there is no $x[0.5]$. It is either $0$ or $1$. Central difference might appear to have a $\frac{1}{2}n$ term but that could be seen as a "phase shift" with respect to the sampling of the signal. The relative temporal difference between two successive samples is still one sampling period.

And this is the interesting bit here. The finite difference is the result of something that happens to two samples. Not one.

So, from the point of view of visualisation, you could show your first order finite difference as "occuring" between samples ($+\frac{1}{2}Ts$, where $Ts$ is the sampling period). The second order finite difference is something that "occurs" between pairs of samples of the first finite difference. The third order finite difference... (and so on).

Notice here that as you are working "up" the finite differences order, the starting point is "sliding" and the total length of the finite difference signal becomes shorter. The first order finite difference is between $x[n],x[n+1]$, the second order finite difference is between $x'[n],x'[n+1]$ which originate from $(x[n], x[n+1])$, $(x[n+1], x[n+2])$, the third order finite difference is between $x''[n], x''[n+1]$ which originate from $(x[0], x[1])$, $(x[1],x[2])$, $(x[3],x[4])$ and so on.

Notice here that throughout this nested evaluation of these finite differences that "I don't know" bit is becoming larger and larger. You cannot evaluate the $y'''[0]$ of the third finite difference if you don't know the values of the first four samples of $x[n]$. But if you consider an "overlapping" segment across a part of the signal that is long enough to evaluate all finite differences, then you would be able to "see" all of them.

Hope this helps.


Deconvolution of a 1D Time Domain Wave Signal Convolved with Series of Rect Signals I have a synthesized signal (the bottom of the following figure), which is the convolution of the input signal (at the top) and the objective function (in the middle). The intention is to retrieve the objective function from the convoluted signal, when the input signal is known. From practical point of view, it seems to be an ill-conditioned problem, but I'm curious to know the expert opinions on how far one can get, and with what SP toolset. Specifically, I would like to know: 1). comparing the first wave packet and second, how to tell whether the first square in the objective function is one with shorter time span but the same amplitude, or with the same time span but smaller-amp; 2). for the 4th square, how to tell it's not two separate short squares rather than a single long one. I've also attached the Matlab codes I used to generate the figures. Thank you very much.

clear; close all;

%% 10MHz incident signal;
fs = 1000e6;
f = 10e6;
wavelength = fs/f;
sig = wavemaker(3.5, f, fs);
figure; subplot(3,1,1);
plot(sig, 'LineWidth', 1); title(strcat('input signal, wavelength =',num2str(wavelength),' data pts'));
axis([0,3300,-2,2])

%% Sparse signal;
N = 3000;                    % N : length of signal
s = zeros(N,1);
k = [50:(50+wavelength*0.1) 500:(500+wavelength*0.6) 1200:(1200+wavelength*1.6) 2200:(2200+wavelength*4.6)];
s(k) = 1;
subplot(3,1,2);
plot(s, 'LineWidth', 1); title('distribution: objective function');
axis([0,3300,0,2])

%% convoluted signal
y = conv(sig,s);
subplot(3,1,3);plot(y, 'LineWidth', 1);hold on;
plot(abs(hilbert(y)), 'LineWidth', 1);
title('convoluted signal between input and districution');
xlim([0,3300])

function x = wavemaker(nCycles, fc, fs)
% function to generate wave packet;
nSample = round(fs / fc * nCycles);
ts      = 1 / fs;
T       = ts * nSample;
t_max   = ts * (nSample-1);
t       = 0: ts: t_max;

x       = sin( 2 * pi * fc .* t);

x = x.*hanning(nSample)';
end

Solving a deconvolution isn't easy even in simulated environment not to mention in practice.
The main trick to solve it is using the proper model / prior for the problem and very good measurements (High SNR).
So basically, for deconvolution we're after:
$$ \hat{\boldsymbol{x}} = \arg \min_{\boldsymbol{x}} \frac{1}{2} {\left\| H \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} $$
Where $ H $ is the matrix form of the convolution operator of the known signal and $ \boldsymbol{y} $ is our measurement samples.
Now, we need to talk about the Convolution form. The convolution has to deal with boundary conditions which might be crucial for good results.
We have few approaches (Basically extrapolation):

Extrapolate with Zeros - Assuming data outside the samples of the signals is zero.
Extrapolate by Nearest Neighbor (Also called Replicate) - The unknown value is extrapolated by the nearest known value.
Extrapolate by Periodic continuation - The data is assumed to be periodic. Hence any missing value is based on that.

The building of $ H $ must match that. Since you used, in your code conv() with no explicit mode it means you basically chose zeros and since your convolution output is full (The default) it means we see the transients as well in the output.
The solution to the above is given by:
$$ \hat{\boldsymbol{x}} = \arg \min_{\boldsymbol{x}} \frac{1}{2} {\left\| H \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} = {\left( {H}^{T} H \right)}^{-1} {H}^{T} y $$
The stability of this solution depends on the Condition Number of $ {H}^{T} H $.
Remark - One could solve this in Frequency Domain as well. Though it will require some touches as in Frequency Domain the model is the periodic model.
Let's have a look on the results:

First we can see the Condition Number is huge!
You may think about the condition number as the amplification of the error. It means even the slightest noise will make things unstable. As can be seen, indeed even a white noise with a standard deviation of 1e-8 caused errors!
In practice, to deal with this instability we use some regularization.
$$ \hat{\boldsymbol{x}} = \arg \min_{\boldsymbol{x}} \frac{1}{2} {\left\| H \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} + \lambda R \left( x \right) $$
Where $ R \left( \cdot \right) $ is the regularization function and $ \lambda $ is the regularization factor which balances between listening to the naïve deconvolution model or  to the regularization model.
The regularization function must be chosen with respect to the prior knowledge we have about the signal of interest.
In your case, something clear about your signal is its piece wise smooth property. Its gradient is very sparse. Hence it is a perfect match to the Total Variation model:
$$ \hat{\boldsymbol{x}} = \arg \min_{\boldsymbol{x}} \frac{1}{2} {\left\| H \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} + \lambda TV \left( x \right) = \arg \min_{\boldsymbol{x}} \frac{1}{2} {\left\| H \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} + \lambda {\left\| D x \right\|}_{1} $$
Where $ D $ is the finite differences operator (Discrete Derivative).
This is a relatively easy task to solve. In my project Total Variation (TV) Regularized Least Squares - Solvers Analysis I implemented and compared few solvers.
For this example I used the ADMM based solver.
Here is the result for the TV regularization:

As can be seen, it totally overcome the (Very Low!) noise from above.
In real world (And higher noise levels) one needs to tweak the $ \lambda $ parameter. You will hardly recover perfect results, but they will be much better than doing the naive approach.
MATLAB Implementation
The full code is available on my StackExchange Signal Processing Q71822 GitHub Repository.
It includes the functions to build the convolution matrix from the samples and solve the TV problem.


Showing linear-phase frequency response when only real part is known I'm studying for an exam, and this particular practice problem is missing the solution:

  Given a causal LTI system where only the real part of the frequency response is known,
  $$\Re\left\{H(j\omega)\right\}= 1 + \cos(\omega)
$$
  show that the system has a linear-phase frequency response.

I know that systems with constant group delay have linear phase, but I can't figure out how to get the group delay without the imaginary part, since group delay is
$$
\tau(\omega)=\Re\left\{\frac{\mathcal F[th(t)]}{\mathcal F[h(t)]}\right\}
$$
and
\begin{align}
\Re\left\{\mathcal F[th(t)]\right\} &= \Re\left\{j\frac{d}{d\omega}H(j\omega)\right\}\\
& = \Re\left\{\frac{d}{d\omega}[jH_{\rm real}(j\omega)+j^2H_{\rm imag}(j\omega)]\right\}\\
& = -\frac{d}{d\omega}H_{\rm imag}(j\omega)
\end{align}

Is my approach right but I made a mistake or can't see some trick, or am I going down the wrong path?

You have two kinds information in the question:

The system is causal
The real part of the Frequency response is given.

Now, (repeating the steps from Wikipedia Entry on Causal Filter)

$h(t)$ is causal, i.e. $h(t)=0, t<0$.
Let $g(t)=\frac{1}{2} (h(t)+h^*(-t))$ which is non-causal, but has hermitian symmetry, hence its Fourier Transform is real.
We know (due to causality) $h(t)=2u(t)g(t)$ where $u(t)$ is the unit-step function
Now, computing Fourier-Transform we get
$$
H(\omega) = 2 \mathcal{F}(u(t))*G(\omega)
$$
where $*$ is the convolution operator.
Now, we can express this as
$$
H(\omega) = \left(\delta(\omega)-\frac{j}{\pi\omega}\right)*G(\omega)=G(\omega)-j\mathcal{H}(G(\omega))
$$
where $\mathcal{H}(G(\omega))$ is the Hilbert Transform of $G(\omega)$. Since we know that $G(\omega)$ is real, $G(\omega)$ equals the real part of $H(\omega)$. Then, $\mathcal{H}(G(\omega))$ is the imaginary part of $H(\omega)$.

So, to get the imaginary part, you need to calculate the Hilbert Transform of $1+\cos(\omega)$ which is given by

$$
\mathcal{H}(1+\cos(\omega))=\sin(\omega).
$$

Accordingly, the phase of your system is given by

$$
\phi(\omega) = \tan^{-1}\left(\frac{-\sin(\omega)}{1+\cos(\omega)}\right)
$$

Now, if we calculate the derivative regarding omega, we need to get a constant for a linear phase filter:

$$\frac{\partial \tan ^{-1}\left(\frac{-\sin (x)}{\cos (x)+1}\right)}{\partial x}=-\frac{\frac{\cos (x)}{\cos (x)+1}+\frac{\sin ^2(x)}{(\cos (x)+1)^2}}{\frac{\sin ^2(x)}{(\cos (x)+1)^2}+1}=-\frac{1}{2}$$

So, this shows that you have indeed a linear-phase system.

Your way with the group delay is also correct, once you have the imaginary part of your $H(\omega)$ obtained from its Hilbert Transform.


Stability of transfer function Now K(s) is obviously the negative feedback loop which is ${H(s) \over 1+H(s)G(s)}$

When I substitute ${H(s) = {1 \over s-2}}$ I get  K(s) = ${ 1 \over s-2 +G(s) }$

For the system to be stable I know that the poles have to be on the negative half of th imaginary plane. Not sure how to find G(s).

Hint:

Since $K(s) = \frac{1}{s-2+G(s)}$, then the roots of $s-2+G(s)$ must lie on the left side of the complex plane. Let's suppose that $G(s)$ is of the form $as + b$. Can you solve it then?


What is 'Normalized frequency in the range [0,1)', à la DTMF & Goertzel algorithm I have a functioning DTMF program. It takes a larger signal and breaks it up into many small components, then analyzes each segment with an FFT to determine if the magnitude of the response for a particular "bin", corresponding to a specific frequency, is present, thereby indicating that the number which is respresented by that frequency has been pressed, straight forward enough, but I want to implement the Goertzel Algorithm of the library jmathstudio and I can't quite figure out how to go about searching for the frequency I want. the documentation says this:

I construct the bins like this, the term '200' is dictated by the fact that my FFT is currently of size 400, and 2000 is determined by the fact that the sample rate is 4000.

int sixNineSeven = (int)(697.0*200/2000);
int sevenSevenZero = (int)(770.0*200/2000);
int eightFiveTwo = (int)(852.0*200/2000);
int nineFourOne = (int)(941.0*200/2000);

int twelveZeroNine = (int)(1209.0*200/2000);
int thirteenThirtySix = (int)(1336.0*200/2000);
int fourteenSeventySeven = (int)(1477.0*200/2000);

I guess the main thing I need to do is this step:

Normalized frequency in the range [0,1)

But does that mean that I can totally ignore the sample frequency?

I suppose I would be doing this Geortzel operation in lieu of an FFT, so should I just search for this normalized frequency in each of the small segments?

But when I do this:

   Complex coeff = su.goertzelFrequencyAnalysis(new Vector(c), (float)(1/697));

   System.out.println("coeff "+coeff);

My output is just a bunch of garbage like this:

   coeff org.JMathStudio.DataStructure.Complex@5e228a02
   coeff org.JMathStudio.DataStructure.Complex@7bd63e39
   coeff org.JMathStudio.DataStructure.Complex@2e8f4fb3
   coeff org.JMathStudio.DataStructure.Complex@42b988a6

The normalized frequency '1' corresponds to the Nyquist frequency, i.e. half the sampling frequency. So the normalized frequencies $w$ are obtained from the frequencies $f$ in Hz by

$$w=\frac{2f}{f_s}$$

where $f_s$ is the sampling frequency.

The frequency (in Hz) corresponding to a certain bin $i$ ($0\le i\le N-1$) of an FFT with length $N$ is

$$f_i=\frac{if_s}{N}\quad i=0,1,\ldots,N-1$$

For real-valued signals, the bins $i=N/2+1$ (assuming $N$ is even) to $i=N-1$ are redundant.


How to calculate the required bandwidth from a spectrum diagram? [FM] For example from this diagram.

There are many definitions of bandwidth; see for example section 2.2 in this book (PDF linked at the end of the page). In your case, it is easy to use the "absolute bandwidth", which is defined as $B=f_{max}-f_{min}$, where $f_{max}$ is the largest frequency in the spectrum, and $f_{min}$ is the smallest.

According to this definition, you'd have $B=(100.028-99.972)\mathrm{\,MHz}=56\mathrm{\,kHz}$.


How to read an image in Scilab In Matlab we use below function to read an image :

var = imread('IMAGE PATH');

But, how to do same thing in Scilab ?

Scilab Image and Video Processing toolbox provides imread() function.

The related How To entry gives example and it's similar to Matlab:

im=imread(filename);

See ATOMS article to learn how to install additional modules in Scilab.


Best method to segment a given image I have one image. I want to segment it by simple approach. Could you suggest to me the best method to segment it? I think one of them is gradient or FCM. But I am not sure. Please see  my given image file here
and segmented image file below

For such a segmentation, you might want to use graphcuts:
http://en.wikipedia.org/wiki/Graph_cuts_in_computer_vision

However lately there have been better algorithms proposed, such as Branch and Mincut:
http://research.microsoft.com/en-us/downloads/994e80f7-007b-4fed-a6e9-db5f4f32ccd1/

For most of the methods, source codes are available and you could try them out of the box. Also check out:
http://vision.csd.uwo.ca/code/


Help understanding FFT Phase: code and results provided I'm trying to figure out how FFT's phase works. Here's my code:

import numpy as np

t = np.linspace(0,120,1200)
acc = lambda t: 10*np.sin(2*np.pi*t)

# signal is a perfect 10 amplitude 1 frequency sin wave
signal = acc(t)

FFT = np.fft.rfft(signal)
freqs = np.fft.rfftfreq(len(signal),1/10.)
real = signal[110]
total = 0
for binnum in range(len(FFT)):
    bin = FFT[binnum]
    amp = (np.sqrt(bin.real**2 + bin.imag**2))/len(freqs)
    angle = np.angle(bin)
    phase = np.arctan2(bin.imag, bin.real)
    freq = freqs[binnum]
    contribution = amp * np.sin(2*np.pi*freq*t[110] + phase)
    total += contribution

    print "Amp: %f, angle: %f, phase: %f, freq: %f" % (amp, angle, phase, freq)

In the print, I see this line:

mp: 9.815303, angle: -1.256637, phase: -1.256637, freq: 1.000000

It thinks the phase of the perfect 1 frequency wave is -1.25, but it should be 0. How do I find the right phase?

Also, I understand if it's because of spectral leakage, but how does the inverse FFT work then? I was under the impression that to reconstruct the signal at time t, you'd just sum up all the contributions from the sinusoidal waves. But it doesn't work, since the output of the above code is:

total: -9.90868489422
real: 0.576119838246

You phase isn't zero because your FFT width isn't an exact multiple of your sinewave's period.