How does a Huffman encoded message gets unified before channel encoding? I aim to understand the signal entire transmission pipeline. Currenlty, I am facing a problem regarding the handling of a Huffman encoded message.

Example:

Message $X$ contains symbols $a(50\%)$, $b(30\%)$ and $c(20\%)$.
Huffman Code with code word length $2$:

\begin{align}
cc &= 11 \\
bc &= 101\\
ac &= 011\\
\vdots &=\vdots\\
aa &= 0000
\end{align}
There are code words with different lengths. When I go further down the signal transmission pipeline i. e. modulation or channel coding, both requisite bit packages of equal length.

Question: Where and how do the Huffman encoded code words get unified in its length?

Suggestion which doesn't work:

I always take 8 bits regardless of the encoded symbol.
$$11111111 = cccccccc$$
$$101101101 (9\text{bits})\implies\text{cutoff 1 bit} \implies 10110110 (8\text{bit}) \implies bcbc?? \\\implies\text{doesn't work}$$

They don't get unified.

Think of the transmitter pipeline (data source, source encoder, channel coder, modulator, etc) as a sequence of independent blocks. Blocks don't assign any particular meaning or order to their input: they regard the input as just a stream of bits.

So, the output of the Huffman encoder can be regarded as a stream of 0s and 1s. The channel coder will subdivide the stream in groups of $k$ bits and convert them to groups of $n$ bits. Again, the channel coder output is a stream of bits that is processed by the modulator.

This approach has several benefits. One is that each block can be studied independently of the others. Also, you can connect, disconnect and replace blocks, and the system will continue to work, because there is no dependency between the blocks.

The main disadvantage is that there is a little bit of extra performance that can be obtained by making the blocks dependant of each other -- see for instance "joint source-channel encoding".

The receiver is more complicated because a block may need to know the boundaries of the code blocks. This is solved by introducing a frame structure to the data.


Hilbert Transform from FFT? For a school project, I'm working on a software-defined radio transmitter intended for the HF amateur radio bands. I'm planning to support SSB transmission with the formula

$$f(t) = m(t) \cos(2\pi f_\text{carrier}t) \pm \hat m(t)\sin(2\pi f_\text {carrier}t)$$

where $f_\text{carrier}$ is the IF carrier and $\hat m$ is the Hilbert transform of $m(t)$. Of course, I will use a discrete form of this equation for my implementation, which means I need to calculate the discrete Hilbert transform of $m(t)$.

I've looked around online but can't seem to find a decent C/C++ implementation of the discrete Hilbert transform that would be suitable for my purposes. However, there are plenty of libraries to calculate FFTs.

From my understanding, a discrete Hilbert transform can be calculated by taking the FFT of the signal and multiplying by j to achieve the 90° shift. It suffers from Gibbs' phenomenon, it seems, and might need a wide bandpass filter.

Can anyone tell me if my understanding is correct (or of a good discrete Hilbert transform function)?

I've looked around online but can't seem to find a decent C/C++ implementation of the discrete Hilbert transform that would be suitable for my purposes.

Can't completely agree to that, see links below.

So, a hilbert filter generates the analytical signal out of a real signal.

In fact, it's but a complex bandpass filter: It cancels out all the negative frequencies.

There's a lot of ways you can implement that (including GNU Radio's hilbert filter, source code, filter generation code), but you can imagine taking a symmetrical half-band low pass filter, and shifting it up for half the sampling rate, so that its passband now covers the positive frequencies.

Luckily, half-band filters are easy to produce and to implement; what is often done is that you take a prototype window, interpolate it by a factor 2 (adding zeros between the window taps), and weigh it, for example with 1/n (n=tap distance from center tap). The result works quite fine:


Wavelet Transform and STFT How wavelet transform is different from STFT.

I'm not able to understand what is resolution in frequency domain means?

In the STFT, you apply windowing and Fourier transform on the signal using sliding patches and then combine the resulting transforms, which will help you eventually end up with a uniform time/frequency representation of the signal.

In the wavelet transform case, you apply a filter bank on the overall signal at once. In this way, you obtain a coarse-to fine resolution pattern on the time/frequency representation.

Both methods result in similar time/frequency representations which can be derived from each other.

The major differences:
(1) STFT is uniform yet CWT is not.
(2) You apply STFT on patches, but you apply CWT on the overall signal.
(3) STFT involves Fourier transforms but CWT only requires an orthogonal filter bank.


A communication system I am solving the following problem as an exercise.
I have a signal
$$z(t) = x(t)\cos(wt) + y(t)\cos(ut)$$
and we assume $x,y$ to be band-limited continuous time signals. In the frequency domain we have that both signals are $0$ if the frequency $\omega \geq \omega_{max}$. We try to reconstruct $x$ by passing $r$ through a bandpass filter $B(\omega)$ that is $1$ if $\omega$ is $\omega = w$ or $\omega = -w$ (it actually is $1$ in an interval centered around $w$ and $-w$ but I didn't know how to write it since I have no information about the width of the interval. The outcome is then multiplied by $\cos(wt)$ and then passed through a lowpass filter $L(\omega)$ which is $2$ if $-\omega_{max}\leq \omega\leq  \omega_{max}$. What do we need to assume on $w,u$ in order for $x$ to be reconstructed?
In the frequency domain, I have that
$$Z(\omega)= \frac12\big(X(\omega+w)+X(\omega-w)+Y(\omega-u)+Y(\omega+u)\big)$$
So picotrially I have 2 copies of $X$ centred at $w$ and $-w$ and two copies of $Y$ centred at $-u$ and $u$ with half the amplitude, right?
Now, the first bandpass filter should give me only the two copies of $X$, since it is centred around $w$, then there is this multiplication with $\cos(wt)$ again. What does this do? Doesn't it simply make yet another two copies of $X$ centred around $w$? Does it give me back the original amplitude getting rid of the $1/2$? Then, how can a lowpass filter allow me to reconstruct the signal when I have $2$ copies of it?

Let me slightly change notation and denote the two modulation frequencies by $\omega_1$ and $\omega_2$.
You need $|\omega_1-\omega_2|>2\omega_{max}$ to be satisfied, otherwise the shifted spectral of $x(t)$ and $y(t)$ will overlap. If this condition is satisfied, you can retrieve $x(t)$ by demodulation (multiplication with $\cos(\omega_1t)$) and (scaled) low pass filtering, as you suggested:
$$\begin{align}2x(t)\cos^2(\omega_1t)=x(t)\big[1+\cos(2\omega_1t)\big]\Big|_{low pass}=x(t)\end{align}\tag{1}$$


Pulse shaping with RRC : Number of taps I am running a BPSK flowgraph in GNU Radio which is based on RRC pulse shape. I noticed that the BER is quite poor when the number of taps is low. However, using the number of taps around $20\times Samples\_Per\_Symbol$ result in performance close to theory. I found the above number through trial and error and I was wondering if there was any general formula for determining the number of taps for RRC pulse shaping.

Regards,

Moses.

See pages 5 and 6 and the plots on the following pages specific to number of samples per symbol in this very helpful reference by Ken Gentile on designing RRC pulse shape filters:

http://www.analog.com/media/en/technical-documentation/application-notes/AN-922.pdf

An example I have previously done shows the consideration of filter length (how many symbols does the filter span) and number of samples per symbol. The top plot shows a 100 tap RRC filter with only 2 samples per symbol and the lower plot it is the same number of taps but 10 samples per symbol.

Observe the following key points that illustrate some of the important considerations involved: With 2 samples per symbol the filter impulse response extends further in time. This results in a greater stop band attenuation for the same filter roll-off factor (Alpha) and complexity (number of taps). The difference is particularly dominant close to the band edge in the area of 700 Hz in the plots.  In the figures shown we see that we get approximately 70 dB of attenuation for 2 samples per symbol but only approximately 35 dB of attenuation for 10 samples per symbol. However the overall frequency of the digital spectrum is most constrained for 2 samples per symbol; this is resolved with subsequent up-sampling or consideration elsewhere in our overall system design. As long as those considerations are properly addressed, according to Nyquist sampling theorem we will achieve the full integrity of our RRC filter design (with only 2 samples per symbol!).
In the second case, because of the higher number of samples per symbol, with the same number of filter taps (complexity) the impulse response of the filter does not extend as far in time, and is therefore truncated earlier. This degrades overall performance as visible here with the inferior filter rejection. However the greater overall digital frequency spectrum is wider.
An efficient design would choose to do two samples per symbol for such a "shaping" filter and then upsample to the bandwidth needed to be compatible with the analog design or other bandwidth considerations, since the upsampling filter design would be relatively simpler. Also to be considered for many applications is a polyphase RRC upsampling filter structure. In this case we would choose an even number of taps and incorporate the upsampling as described further in this response: How to implement Polyphase filter?

Keep in mind, as explained further in this post: Why root raised cosine filter can eliminate intersymbol interference (ISI) ? the entire motivation to do RRC filtering is spectral efficiency: hence regulatory requirements (spectral mask) will often dictate the rejection ultimately needed in the RRC filter (as long as subsequent analog stages, in particular the power amplifier, maintain sufficient linearity to not undo the shaping we so carefully applied...resulting in "spectral regrowth!).


Meaning: channel unknown and known at the transmiiter What does "channel unknown" and "known at the transmitter" mean?

In books about signal processing, it has always mentioned it if capacity or transmission (estimation) is discussed. I can't understand the difference.

If I simulate the transmission, the channel is given to me (the channel is known at the Transmitter). I can model it using the standard function in Matlab or as a random function.  In estimation algorithms, like ZF, it is also used as a channel ( inverse of channel matrix) for an estimate. What is the difference? How can channel be unknown?

edit 1.

It is unclear what difference is between CSI unknown and known. If I  simulate transmission system, I need a channel. If H is given as a matrix or via rayleigh-function, it means CSI is known, isn't'?  The Capacity is computed using a channel matrix (in its equation).  I cant understand then why there are two difference equation: CSI known and unknown to the transmitter.

You need to distinguish between transmitter (TX) and receiver (RX) here. For the receiver it is generally quite easy to gain knowledge of the channel: have the transmitter send a few known symbols, then the receiver can estimate the channel coefficients (e.g., via Least Squares or MMSE). The receiver needs this channel estimate for equalization. For example, the ZF receive strategy you mentioned essentially means you divide out the channel coefficient, i.e., if you receive $y = h \cdot s + w$ where $h$ is the channel, $s$ your unkown information symbol and $w$ the noise, the simplest you can do to estimate $s$ is divide out $h$ (use $y/h$ as an estimate), for which you need to know it.

The question you asked about is knowing the channel at the transmitter. This is much harder since in general, the receiver would need to send its estimated channel coefficient back to the transmitter to know it. This requires extra communication resources. And time, which is a problem since the channel information  might be outdated by the time it reaches the transmitter, if your channel is time-varying. Another option is to assume the channel between TX and RX are "symmetric", in communication parlance referred to as reciprocity. Then, the TX can estimate its own received channel for data it got from the RX and assume that the same channel coefficient applies for the reverse link. For this to work, in general we need to assume that TX and RX work on the same frequency (e.g., when they use time division duplexing).

The real question here is what would it help the TX to know the channel? Does that improve capacity? As I'm sure you read it does help in the MIMO case as this allows the TX to send data in the "proper directions", i.e., excite the strong eigenmodes of the channel and thus get more power across. There are other cases where it helps, e.g., when dealing with multiple users, interference, and other things.


Wiener-Hopf equation The Wiener-Hopf is given by

$$\mathbf W_\text{opt} = \mathbf R^{-1} \mathbf P$$

What if the matrix $\mathbf R$ is almost singular?
When $\mathbf R$ is singular, the inverse of $\mathbf R$ does not exist,
what modification can I do to make it usable?

My idea is move the $\mathbf R$ to the left side of the equation, is that right?

I guess, you refer to the Wiener-Hopf equation in discrete time for FIR filters.. If this is a homework-task, then I guess most probably you did some mistake in calculation of $R$ (e.g. if you have white noise in the signal, you wont have a singular $R$).

Anyway, there are some methods to solve this:

You dont necessarily need the inverse of $R$. You need a solution to the problem $RW_{opt}=P$, which leads to $W_{opt}=R^{-1}P$ in case $R$ is regular. However, if $R$ is singular, it may still happen, that $P$ is in the image of $R$. Then, still you can find $W_{opt}$ as the solution to the equation $RW_{opt}=P$.
(this is just an approximate solution): If you have a singular $R$, you can replace the equation by $W_{opt}=(R+\alpha I)^{-1}P$, i.e. add some regularization factor to the matrix. This will make the system well-defined, but it will not correspond to the exact solution. Note, that in case $R$ is singular and $P$ is not in its image, then there is no exact solution.
You can go for the pseudo-inverse of $R$. BUt again, this will not correspond to an exact solution, as it does not exist.

All in all, if your $R$ is singular, something is really bad in your model I believe.


Fast overlap save method If x[n]= (some values), h[n]=(some values).  Now suppose length of h[n] is a multiple of two. Being specific say it has length of four.. Then M=4, Then what should be value of L and value of N to satisfy N=M+L-1..

E.g. which i tried solving is.

x[n]={7,6,4,5,2,4,5,2,3}
h[n]={1,2,3,1}

Taking N=8 i got L=5.

x1[n]={0,0,0,7,6,4,5,2}
x2[n]={4,5,2,4,5,2,3,0} (extra zero to make length=N)

Now while solving I found to lose some of values using save method.

for overlap-save, (sometimes called "overlap-scrap") when you (circularly) convolve h with x1. you scrap the first 3 output values of y1 and keep the last 5 values.  when you (circularly) convolve h with x2, again you scrap the first 3 values of y2 and keep the last 5 values.

with

x[n]={7,6,4,5,2,4,5,2,3}
h[n]={1,2,3,1}

it should be

x1[n]={0,0,0,7,6,4,5,2}
x2[n]={4,5,2,4,5,2,3,0}

remember that in convolution, h[n] gets turned around to be h[-n] (h[0] remains in the same position but the other values are extended to the left).  only when h[0] is aligned with x1[3] will you get a legitimate value for y1[3].  and notice, because you (necessarily) delayed x[n] by 3 samples with that initial zero-padding, the output will also be delayed by the same amount.


Cepstrum, peak picking I try to obtain fundamental frequency from cepstrum. I decide to use peak picking algorithm, which pick maximum value of the Cepstrum (c) with sampling rate fs.

function [f0] = spPitchCepstrum(c, fs)
 % search for maximum  between 2ms (=500Hz) and 20ms (=50Hz)
 ms2=floor(fs*0.002); % 2ms
 ms20=floor(fs*0.02); % 20ms
 [maxi,idx]=max(abs(c(ms2:ms20)));
 f0 = fs/(ms2+idx-1);
end

It works but I can not understand why f0 is calculated like:

f0 = fs/(ms2+idx-1);

where in denominator is value consist indexes.

Can you somebody explain me why it is calculate in this way?

Cepstrum argument is called quefreency, which is in fact a time domain. So for example if you are looking for the fundamental frequency then you are searching for a peak in a specific range. In your case that's $[0.002; 0.2] \;sec$, which corresponds to $[50; 500] \; Hz$., knowing that $f=\dfrac{1}{t}$.

So when you are searching for fundamental frequency, you can find peak in quefrency domain and then convert it to corresponding frequency. You can also use a different approach. We know that our cepstrum is symmetrical around center (similarly to Fourier Transform), and also we have our sampling frequency $f_s$. Let's say that you have $N$ cepstrum bins, therefore distance between each of them is: $\dfrac{f_s}{N}$. That's why you can find the corresponding fundamental frequency $f_0$ by dividing sampling frequency by the index of a peak.

For example our peak is at index number $10$, and sampling frequency is $1000 Hz$. Corresponding $f_0$ is $100 Hz$. If peak is at bin $50$ (more to the right, meaning lower frequency - remember quefrency is like time), then $f_0 = 20$ Hz. The subraction of 1 in denominator is cause by fact that MATLAB is indexing from 1 not from 0, like most of the programming languages.


The use of LS estimation with Vector OFDM I work in Vector OFDM, but I had an issue in channel estimation in that scheme. For example, I try to use LS channel estimation but I can't get the performance back. However, in conventional OFDM, I can use that estimation straightforward.

First, let me present the difference of VOFDM compared with OFDM, the V-OFDM scheme is presented as below figure:

As you see, the only difference that in V-OFDM divides the data bock into vectors and then process them by vectors.

My question is how to use channel estimation based on pilot insertion in this case ?? any channel estimator such as for example LS. I tried different inserting the pilots in similar way of that in OFDM, but I can't get the performance back !!

Thank you

Sorry for late reply .. I was little bit busy.

you have mistakes in your code. Although mythology is right, you have mistakes in some parameters.

Check this paper "Low-Complexity Equalization of Orthogonal Signal-Division Multiplexing in Double-Selective Channels" Then,  Modify your code following part II in that paper. It's with details there, if you got any issue, you'd let me know.
Regarding the channel estimation, you follow the same way of conventional OFDM scheme except for V-OFDM you use Vectors pilots instead of symbol pilots.  I means in that mentioned paper, you will have vector of size $M$ used as pilots.

Good luck


Fourier transform relationship I am having trouble understanding the relationship between a frequency function and it's Inverse Fourier transform.

The Frequency function is

$$\frac{1+0.8(e^{-j 2\pi f}+e^{j 2\pi f})+0.64}{1+1.4\cos(2\pi f)+0.49}$$

And the inverse Forier transform is
$$\frac{(-0.7)^{|n|}}{1-0.49}+\frac{0.64(-0.7)^{|n|}}{1-0.49}+\frac{0.8(-0.7)^{|n-1|}}{1-0.49}+\frac{0.8(-0.7)^{|n+1|}}{1-0.49}$$

Consider the signal $$s[n]=\frac{a^{|n|}}{1-a^2},\, |a|<1$$

Assuming $$x[n]=a^n u[n]\stackrel{\text{DTFT}}\longleftrightarrow\frac{1}{1-a e^{-j\omega}},$$
 one can show that

$$s[n]=x[n]*x[-n]$$

The proof is simple. For example, for $n>0$ we have
$$\sum_{k=n}^{\infty}x[k]x[k-n]=\sum_{k=n}^{\infty}a^ka^{k-n}=\sum_{k=n}^{\infty}a^{2k-n}=a^n(1+a^2+a^4+\cdots)=\frac{a^n}{1-a^2}$$

Hence, the DTFT of the signal $s[n]$ is given by
$$S(e^{j\omega})=X(e^{j\omega})X(e^{-j\omega})=\frac{1}{1-a e^{-j\omega}}\frac{1}{1-a e^{j\omega}}=\frac{1}{1-2a\cos(\omega)+a^2}$$
In summary,
$$\boxed{\frac{a^{|n|}}{1-a^2}\stackrel{\text{DTFT}}\longleftrightarrow\frac{1}{1-2a\cos(\omega)+a^2},\, |a|<1}$$

Now consider

$$G(f)=\frac{1}{1+1.4\cos(2\pi f)+0.49}=\frac{1}{(1+0.7e^{-j2\pi f})(1+0.7e^{j2\pi f})}$$

whose IDTFT can be calculated using the above formula for $a=-0.7$ as

$$g[n]=\frac{(-0.7)^{|n|}}{1-(-0.7)^2}$$

Since
$$\frac{1+0.8(e^{-j 2\pi f}+e^{j 2\pi f})+0.64}{1+1.4\cos(2\pi f)+0.49}=G(f)+0.8e^{-j 2\pi f}G(f)+0.8e^{j 2\pi f}G(f)+0.64G(f)$$
The IDTFT is
$$g[n]+0.8g[n-1]+0.8g[n+1]+0.64g[n]$$
But obviously, the variable $\tau$ is used instead of $n$ in the give equation.


Inverse fourier transform of Hermitian function, getting an imaginary part in the cartoon below

It shows that if we take the inverse Fourier transform of a Hermitian function, real part even and imaginary part is odd we should get a purely real function in the time domain.

I tried to replicate this by taking a frequency response I have, zero padding it making it even for the real and odd in the imaginary. Shown below

Although my ifft in Matlab I still have a significant imaginary component, the real is shown left and imaginary on the right.

I also wrote my own script for DFT/IDFT and tested it on a sine function to make sure I get the same as Matlab. Maybe I have missed something?

A (near) sinusoid in one domain (zoom in off center, and you might see it) can indicate a shift away from symmetry in the other domain.  So your initial imaginary signal may be 1 or 1/2 samples off from exact circular symmetry around element 0 (maybe around element 1 in MatLab or Fortran weird array indexing)


How to calculate the impulse response from the transmission function of an optical device? Let $f(x,y)$ be the electric field in complex representation processed by an optical device described by a transmission function $t(x,y)$ such as:
$$
f'(x,y)=f(x,y)t(x,y)
$$

Where $f'$ is the output electric field.

How can I calculate the impulse response of the device? I don't understand how to determine a function $h(x,y)$ such as:
$$
f'(x,y) = f(x,y)*h(x,y)
$$

More specifically I'm working with a lens transmission function in the paraxial approximation:
$$
t(x,y)= P(x,y)\exp\left[-\frac{jk}{2f}(x^2+y^2)\right]
$$
where $P(x,y)$ is the pupil function of the lens (usually a $circ(x,y)$ function).

P.S.
I wonder if there is a general solution to the mathematical problem of finding a function $h$ for a given $t$ which satisfy $ft=f*h$ for a generic function $f$.

I've posted the same question here: https://physics.stackexchange.com/questions/304444/ but I haven't received any answer.

If the function $t(x,y)$ is an optical transfer function (OTF) then your goal is not to find a function $h(x,y)$ such that the convolution with your input $f(x,y)$ equals the multiplication with $t(x,y)$. Convolution with the impulse response (point spread function (PSD) in optics) happens in a different domain. The OTF is the Fourier transform of the PSD, and the equality

$$f'(x,y)=f(x,y)t(x,y)$$

corresponds to the equality

$$\mathcal{F}^{-1}\{f'(x,y)\}=\mathcal{F}^{-1}\{f(x,y)\}\star h(u,v)$$

where $h(u,v)$ is the impulse response (PSD) given by

$$h(u,v)=\mathcal{F}^{-1}\{t(x,y)\}$$


What techniques are available for correcting constellation rotation due to phase ambiguity? In order for a receiver to recover a baseband signal, it needs a carrier signal equivalent to that used by the transmitter for upconverting the original baseband signal. Unless pilot tones are used, a receiver doesn't have the original carrier available and has to extract it from the signal itself, with a technique such as a phase locked loop.

However, the phase of the carrier signal to be recovered is often ambiguous to the receiver: Unless there is some asymmetry in the transmission (rotationally asymmetric constellation, non-uniform distribution of symbols), the recovered carrier may end up at one of several phase offsets.

A typical 4-way symmetric constellation like QPSK causes a phase offset of nπ/2 radians where n may be 0, 1, 2 or 3. This phase offset then rotates the received symbols around the origin in one of four different ways:

What techniques are available for decoding the symbols regardless of this rotation? I can think of seven, but all have drawbacks:

Transmit a pilot tone (the original carrier, or a signal derived from it), eliminating the phase ambiguity. Wastes bandwidth and power.
Use an asymmetric constellation which aids the carrier recovery algorithm, eliminating the phase ambiguity. More complex to decode, less bandwidth/power efficient.
Use a symbol mapping where the symbols reside on a rotationally symmetric path around the origin, and encode the data as the distance along that path between successive symbols. Otherwise seems like an elegant solution, but any single incorrectly received symbol causes two errors.

Every once in a while, transmit some unique sequence ("mid-amble") which allows the receiver to easily de-rotate the constellation. Slower aquisition, decreases bitrate.
If the bit error rate is high, rotate the constellation after X symbols. Slower aquisition.
Decode the same constellation in all possible ways at the same time, and pick the output with the lowest bit error rate. Wasteful of processing resources.
Use a symbol alphabet designed in such a way that it works with every possible phase offset. Large overhead, which reduces the bitrate.

While some of these approaches are less practical than others, none are completely satisfactory. How is this issue typically handled? Can the channel coding be designed to be rotation-agnostic without sacrificing bitrate?

In my experience from the space industry, both near-earth and deep-space communication system resolve phase ambiguity by means either unique word detection or different flavors of differential encoding. In unique word detection technique, an "Attached Sync Marker (ASM)"(a.k.a a unique word), a preamble of some sort (usually the hex pattern 0x1ACFFC1D) is attached at the beginning of every frame (known as channel access data unit (CADU)). Now, the receiver rotates the received symbols until the ASM is detected. The diagram below from CCSDS Radio Frequency and Modulation Systems--Part 1: Earth Stations and Spacecraft, page 2.4.11 gives a summary of scenarios for (offset) QPSK where the techniques are more appropriate.

Regards,

Moses.


Extra frequencies after applying low pass filter I used Python/scipy to create an input data (sinewave with frequencies {600,1200,1800,2400,...,9000}) and sampling at 10000Hz. A lowpass filter is created using butter(10, 4000/(10000/2), 10000, btype='low'). After filtering the filter using lfilter function, I plotted the FFT of the input data the filtered data, I found some extra frequencies in the FFT graphs. I think it should show the frequencies at {600, 1200,...,4800}. Where are these extra frequencies from? and how to remove them?

The issue is that you're violating Nyquist which requires that the sampling rate be at least twice the rate of the highest frequency in the signal. In your case all of the sines that you're generating above 5000 Hz are folding back. For example, a 5100 Hz tone (5000+100) would fold back to 4900 Hz (5000-100) and so on.

The 1 kHz tone you highlighted is the 9 kHz tone you generated, but folded back. If you look at your graph more closely you'll notice that this tone is in the unfiltered signal and is in no way related to the filtering.


dft of sampled sine using python I'm trying to write a python script to perform a 100-point DFT over a finite length sample of a sinewave at 1/8 the sampling frequency.  I'm curious why my DFT magnitude plot has two spikes in it when i plot it.

(I'm wondering if k=0 corresponds to $\omega=0$, k=100 corresponds to $\omega=2\pi$, and k=50 corresponds to $\omega=\pi$... and i'm just seeing an alias imagine above k=50.. or if i'm doing something wrong with DFT calculation..)

# Python 3.5 script
import sys
import numpy as np
import matplotlib.pyplot as plt

pi = np.pi

def stem_xy(x, y, title, filename):
    print("=> STEM_XY")
    fig = plt.figure()
    plt.title(title)
    plt.stem(x, y)
    print("creating file: ", filename)
    fig.savefig(filename)
    plt.close()

def plot_xy(x, y, title, filename):
    print("=> PLOT_XY")
    fig = plt.figure()
    plt.title(title)
    plt.plot(x,y)
    print("creating file: ", filename)
    fig.savefig(filename)
    plt.close()

def dsp_dft_magnitude(DIN, N, filename):
    print("=> DSP_DFT_MAGNITUDE")
    #N: Number of bins
    w_0 = 2*np.pi / N #Fundamental Frequency (Radians)
    X   = DIN
    XK  = np.zeros(N, dtype=np.complex)
    XM  = np.zeros(N)

    print("N  = ", N)
    print("W0 = ", w_0)

    for k in range(0,N):
        XK[k]=0
        for n in range(0,len(X)):
            XK[k] = XK[k] + X[n]*np.exp(-1j*k*w_0*n)
        XK[k] = (1/N)*XK[k]
        XM[k] = np.sqrt((XK[k].real)**2 + (XK[k].imag)**2)

    k = range(0,len(XM))
    plot_xy(k, XM, "DFT", filename)

    return XM

def run():
    print("=> RUN")

    sampleF   = 1E3
    sampleN   = 100
    Fsin      = sampleF/8

    sampleT   = 1 / sampleF
    t         = np.linspace(0.0, sampleN*sampleT, sampleN)
    y         = np.sin(2*pi*Fsin*t)

    stem_xy(t, y, "samples", "plot_test0.png")
    dsp_dft_magnitude(y, 100, "plot_test1.png")
    sys.exit()

run()

Found the answer in numpy documents for fft:

# python to perform dft
# from import numpy.fft import *

A = fft(a, n)

  A[0] contains the zero-frequency term (the sum of the signal), which is always purely real for real inputs.

  A[1:n/2] contains the positive-frequency terms

  A[n/2+1:] contains the negative-frequency terms, in order of decreasingly negative frequency.

  For an even number of input points, A[n/2] represents both positive and negative Nyquist frequency, and is also purely real for real input.

  For an odd number of input points, A[(n-1)/2] contains the largest positive frequency, while A[(n+1)/2] contains the largest negative frequency.

There are many ways to define the DFT, varying in the sign of the exponent, normalization, etc. In this implementation, the DFT is defined as

$$
A_k = \sum_{m=0}^{n-1} a_m exp \left\{ -2 \pi i \frac{mk}{n} \right\},\ \ \ \ k=0,...,n-1
$$

The DFT is in general defined for complex inputs and outputs, and a single-frequency component at linear frequency  is represented by a complex exponential , where  is the sampling interval.

The values in the result follow so-called “standard” order:

Correction to code above:

# python DFT Sin Test
# Bill Moore
# 4/11/2019

import sys
import numpy as np
import matplotlib.pyplot as plt

def dsp_mag_and_phase(DIN, N, filename):
    """
    Cleanup and plots DFT and FFT results

    N        : Number of frequency bins
    filename : png image file to save dft magnitute plot
    """

    print("=> DSP_FREQPLOT")

    XK = DIN
    K  = range(0, N)

    # Swap Positive and Negative Frequencies of DFT/FFT
    #  (to get standard spectrum of -pi to pi)

    # Assign a positive or negative frequency to each point
    KK = []
    for k in K:
        if k>=0 and k <= int(N/2):
            KK.append(k)
        else:
            KK.append(k - int(N))

    #print("KK:", KK)
    #print("len(KK)=", len(KK))
    #sys.exit()

    AM = np.zeros(N)
    AP = np.zeros(N)
    AM = np.zeros(N)
    AP = np.zeros(N)

    #Positive Frequences
    AX_pos = XK[0:int(N/2)+1]
    AK_pos = KK[0:int(N/2)+1]

    #Negative Frequences
    AK_neg = KK[int(N/2)+1:]
    AX_neg = XK[int(N/2)+1:]

    # Concatenate negative then positive
    AX = np.concatenate((AX_neg, AX_pos))
    AK = np.concatenate((AK_neg, AK_pos))

    # Convert Complex Frequency into Magnitude and Phase
    for k in range(0, len(AK)):
        AM[k] = np.sqrt((AX[k].real)**2 + (AX[k].imag)**2)
        AP[k] = np.arctan2(AX[k].imag, AX[k].real)

    return (AK, AM, AP)

def dsp_dft(DIN, N, filename):
    """
    plots N-point DFT

    N        : Number of frequency bins
    filename : png image file to save dft magnitute plot
    """

    print("=> DSP_DFT_MAGNITUDE")

    DFT = np.zeros(N, dtype=np.complex)
    K   = range(0, N)
    w0  = 2 * np.pi / N #Fundamental Frequency (Radians)

    # DFT
    for k in K:
        DFT[k]=0
        for n in range(0,len(DIN)):
            DFT[k] = DFT[k] + DIN[n]*np.exp(-1j*k*w0*n)
        DFT[k] = (1/N)*DFT[k]

    # Covert Complex to Magnitude and Phase
    (AK, AM, AP) = dsp_mag_and_phase(DFT, N, filename)

    plot_mag_and_phase(AK, AM, AP, "DFT", filename + ".png")

    return (AK, AM, AP)

def dsp_fft(DIN, N, filename):
    """
    plots N-point FFT
    (FFT is just a computationally faster DFT)

    N        : Number of frequency bins
    filename : png image file to save dft magnitute plot
    """

    print("=> DSP_FFT_MAGNITUDE")

    # FFT Points equal to length of DIN (default for fft function)
    if N == 0:
        N = len(DIN)

    FFT = np.fft.fft(DIN, N)

    (AK, AM, AP) = dsp_mag_and_phase(FFT, N, filename)

    plot_mag_and_phase(AK, AM, AP, "FFT", filename + ".png")

    return (AK, AM, AP)

def plot_stem(x, y, title, filename):
    print("=> PLOT_STEM")
    fig = plt.figure()
    plt.title(title)
    plt.stem(x, y)
    print("creating file: ", filename)
    fig.savefig(filename)
    plt.close()

def plot_mag_and_phase(k, mag, phase, title, filename):
    print("=> plot_mag_and_phase")

    fig = plt.figure()
    plt.clf()

    plt.subplot("211")
    plt.title(title)
    plt.ylabel("Magnitude")
    plt.plot(k,mag)

    plt.subplot("212")
    plt.ylabel("Phase")
    plt.plot(k,phase)

    print("creating file: ", filename)
    fig.savefig(filename)
    plt.close()

def dsp_range(sampleN, sampleF):
    """
    return  : time points for sampling at a given frequency

    sampleN : Number of Samples
    sampleF : Sample Frequency in Hz
    """
    sampleT       = 1 / sampleF
    time_range    = np.linspace(0.0, sampleN*sampleT, sampleN)
    return time_range;

def run():
    print("=> RUN")

    sampleF   = 1E3
    sampleN   = 100
    Fsin      = sampleF/4

    # sample function
    tvec       = dsp_range(sampleN, sampleF)
    y1         = np.sin(2*np.pi*Fsin*tvec)
    plot_stem(tvec, y1, "samples", "plot_test0.png")

    dsp_dft(y1, 100, "plot_test1")
    dsp_fft(y1, 0,   "plot_test2")

    sys.exit()

run()


differential phase tracking loop for the receiver For a DBPSK modulation, how can we track the differential phase? I thought about using a phase tracking loop but I did not find articles for this.

There is no need to track the phase difference at all.  What is needed is to compare the "average" RF phase during the current bit interval to the "average" RF phase during the previous bit interval.  If the two phases are the same, or near enough the same, the receiver decides that the phase did not change during the transition from the previous bit interval to the current bit interval and so the differential encoding at the transmitter is telling us that a $0$ data bit is being transmitted.  If the two phases are vastly different, then we are being told that a $1$ data bit is being transmitted.

How is all this done? Well, use the part of a QPSK receiver that separates the RF signal into I and Q signals, and then match-filter (or otherwise process) the I and Q signals to produce complex-valued outputs $Z_1 = (X_1,Y_1)$ and $Z_2 =(X_2,Y_2)$
where $X$ is
the output of the I branch and $Y$ the output of the Q branch of the receiver. Note that the $X$'s and $Y$'s are individual real numbers: they are the sample values obtained from the I and Q branch filters/processors at the chosen sampling instants, not continuous-time signals or vectors of sample values that DSP processors or computers use to represent continuous-time signals.

The  hard-decision DBPSK decision device considers the question:

  Is the new symbol $Z_2$ closer to the old symbol $Z_1$ or closer to the
  negative $-Z_1$ of the old symbol?

and thus compares

$$(X_2-X_1)^2 + (Y_2-Y_1)^2 \gtrless (X_2+X_1)^2 + (Y_2+Y_1)^2$$

which can be simplified to a sign comparison on $\langle Z_1,Z_2\rangle = X_1X_2+Y_1Y_2$.
Note that this is essentially asking

  Are the two vectors $Z_1$ and $Z_2$
  are pointing in roughly the same direction (in which case the inner product
  or dot product is positive) or in roughly opposite direction (in which case
  the dot product is negative)?

A third viewpoint thinks of $Z_1$ and $Z_2$ as
complex numbers and asks

  Is $\text{Re}(Z_1Z_2^*) = X_1X_2+Y_1Y_2$ positive or negative?

But, regardless of the viewpoint, the end result is the same.

Note that there is no explicit phase comparison in all this: those boxes marked $\arctan\left(\frac yx\right)$ on block diagrams of communications receivers are more for conceptual understanding than for implementation guidance.

Some notes to read before downvoting this answer.

It is assumed that the receiver has a good estimate of the received carrier frequency but not of the carrier phase.  Frequency acquisition can be done by squaring the RF signal (to eliminate the phase modulation), high-pass filtering (to remove the DC in the squared signal which has DC plus an unmodulated double-frequency signal in it), and frequency-halving, to recover the carrier frequency. In fact, the output of the frequency halver can be used as the local oscillator in the receiver front end.  It is not necessary to track the RF phase.  All that we are interested in is whether the phase in the current interval was the same as the phase in the previous interval or was it different? Where exactly $Z_1$ and $Z_2$ individually lie in the complex plane is not of importance; their relative position is.  This is different from phase-coherent demodulation of plain vanilla BPSK in which there is a single (real-number) output that must be compared to $0$ to determine the transmitted bit and woe betide the designer who ignores the possibility of phase ambiguity caused by the frequency halver.  Note also that we are assuming that the bit interval timing is known with a fair degree of accuracy so that the DBPSK receiver knows when it is the correct time to sample the I and Q branch outputs to get the $Z$'s.


Estimate the Discrete Fourier Transform / Series of a Signal with Missing Samples Assuming we have a discrete signal $ { \left\{ x \left[ n \right] \right\}}_{n = 1}^{N} $.
Which has a Discrete Fourier Transform / Series.

Now, assume I'd like to estimate its Discrete Fourier Series coefficient yet some samples of $ x \left[ n \right] $ are missing (The indices are known).

How could that be done efficiently without computing the Pseudo Inverse of the adapting Fourier Series matrix?

Given $ \left\{ x \left[ n \right] \right\}_{n \in M} $ where $ M $ is the set of indices given for the samples of $ x \left[ n \right] $.

The trivial solution (Which it would be great to have a faster more efficient solution is what I'm looking for) would be:

$$ \arg \min_{y} \frac{1}{2} \left\| \hat{F}^{T} y - x \right\|_{2}^{2} $$

Where $ \hat{F} $ is formed by subset of columns of the DFT Matrix $ F $ matching the given indices of the samples, $ x $ is the vector of the given samples and $ y $ is the vector of the estimated DFT of the full data of $ x \left[ n \right] $.

The solution is then given by the Pseudo Inverse (Least Squares Solution):

$$ y = { ( \hat{F} \hat{F}^{T} ) }^{-1} \hat{F} x $$

In practice, the matrix will be very poorly conditioned hence solution must be generated using the LS Solution using the SVD.

A sample code is shared on GitHub Repository.

Result of the code:


The effects of frequency resolution in FFT on the frequency spectrum This is a section for the parameter study of my thesis (basically: why and how I did what I did in the experiments). Currently, I have to choose the fit bandwidth and FFT lines combination for the experiment (for the parameter study, it is just giving an impulse for the generator, or hitting the metal specimen with a hammer). Investigations (actual measurement) say that with frequency resolution under 0.1 Hz, the spectrum obtained would be very squishy. For example, right below is 12800 FFT lines (or 12800 FFT bins) with different bandwidths (0.5, 1, 2, 5, 10 and 20 kHz) and different frequency resolutions (39.1, 78.1, 156.3, 390.6, 781.3 and 1562.5 mHz respectively)

The spectrum with 39.1 and 78.1 mHz frequency resolution (blue and orange) are very "squishy"...
Question 1: Why do we have such squishy spectrum?
For reference: this is non-squishy spectrum

Question 2: At low-frequency range (<50 Hz), smaller FFT bin shows peak better (more data point on the same range of frequency, thus easier to spot the peaks). But why at high-frequency, smaller FFT bins lead to lower peaks (at 300 and 330 Hz)?
Question 3: Is there a "critical point" for frequency resolution? Or the trade-off between frequency resolution vs spectral/time resolution? If yes, is there any source you can point towards for me, or is this more like a rule of thumb?
Thank you for your support
(This question is also posted in Electrical engineering, and per the first comment, it is put here)

What may explain why your spectrum is noisy is that you are computing it using a single burst of data. You will have to smooth it by averaging successive spectrums applied on your measurements. The number of samples that shall be used for the averaging have to be tuned depending on the speed of variation of your phenomena. After having done this, the other questions can be addressed more easily.
Best regards,
Mourad


Cannot understand the attenuation in signal as the audio bit depth of adc changes I am working on Active Noise Cancellation Project. I came across this strange problem, where the signal is getting attenuated as the audio bit depth changes and I cannot understand it!
My circutit :I connected Audio Output( Head phone terminal) of PC to Audio IN of NI myRIO 1900 device using AUX cable (two way 3.5mm jack cable). I have a pre-recorded signal(.wav) of sinusoidal waveform with Amplitude 1, Frequency 2kHz, 40000 Sample Rate, 80 Number of Samples, Bit depth( bits per sample) = 16.
Now, I just played .wav file using windows media player and tried to record the waveform on NI myRIO module. Surprisingly, I am getting an attenuation of the signal corresponding to the volume reduction in PC(not linear!!). Only I could see the waveform close to amplitude 1, when I increased the volume of the system to 100%. Please explain this !!

The ADC of NI myRIO has its Resolution as 12 bits but the audio Output(.wav file) is 16 bit data. I am expecting it does not have to matter except in quantization errors, as at the PC Audio Out, the audio data should be converted to analog and and again at the 'NI myRIO Audio In', this analog signal will be sampled! (please correct me if I am wrong !)
Another question is, I have tried the same with sinusoidal sound of Amplitude 2.5 and all other specifications are same as given above, now the audio signal is clipped off at 1V above and below -1V.

 ADC(NI myRIO) nominal range is 2.499V to -2.5V! Is it bacause of the bit depth conversion?
Please explain these things !TQ

From NI myRIO-1900 User Guide and Specifications:

  The audio inputs are left and right stereo line-level inputs with a
  ±2.5 V full-scale range

  • For Audio In/Out,
  LSB Weight = 5 V ÷ 212 = 1.221 mV
  Maximum Positive Reading = +2047 * 1.221 mV = 2.499 V
  Maximum Negative Reading = -2048 * 1.221 mV = -2.500 V

Your computer's head phone output probably has another constant of proportionality between voltage and numerical value. Perhaps you can find manufacturer's data, but it is probably safer to simply measure with the help of your NI device.


Fitting a polynomial to a ridge in an array I have an array of numerical data. Say,
%% setup an array of gaussians where the mean is quadratically swept
p = [0.15,-0.3,1];
mu = polyval(p, 1:8);
sig = 2;
tmp = (1/(sig*sqrt(2*pi)))*exp(-(1/2)*(((1:8)'-mu)/sig).^2);
%% add interference
tmp(6,5) = 1.2/(sig*sqrt(2*pi));
%% add noise
tmp = tmp + 1.5e-2*randn(size(tmp));

Assuming that the gaussians are oversampled (sample maxima are close enough to true maxima).
I would like to estimate p from the output array, p_hat.
I have been looking at it for some time, and I cannot twist my head around how to relate this to stuff that I know.

The most straight forward solution would be to take the maximum of each column, do a 1-d fitting using ordinary least squares, solve the outlier somehow. But then the column containing the outlier is not counted at all. I was hoping to be smarter than that, to improve the fit by depending on the secondary (true) peak of the interfered column.
On the face of it, this sounds related to ordinary least squares. But there are 2 independent variables (x-dim and y-dim). And the data is not a sparse point-cloud of "binary" data, but rather a uniform sampled array with amplitude.
The thought to traversing the array left to right (making assumptions of continuity and start or end-point) in a linear-programming kind of way occured to me, but then I would be fitting a continous line with constraints, rather than a parametric polynoma.

I prefer to not use toolboxes. Thanks for any ideas!

Just a couple of thoughts:

That sounds like the classical application for a Kalman filter: you have one gaussian-overlaid signal that "develops" over the columns
If you have a signal model for your signal of interest, a parametric estimator might indeed be the best choice here.
You could write a likelihood function for each observation column $\mathbf c$, like $f_{\mu|\mathbf c}(\mu)=\ldots$ that describes how likely a given $\mu$ is given the observation in each column.
Then, you could state $\hat \mu = \arg\max\limits_\mu f_{\mu|\mathbf C}(\mu)$.
If you can assume your columns to be independent: $f_{\mu|\mathbf C}=\prod\limits_{\mathbf c \in \mathbf C}f_{\mu|\mathbf c}$, so that
\begin{align}\hat \mu &= \arg\max\limits_\mu f_{\mu|\mathbf C}(\mu)\\
&= \arg\max\limits_\mu\;\log\left(f_{\mu|\mathbf C}\right)\\
&= \arg\max\limits_\mu\;\log \left(\prod\limits_{\mathbf c \in \mathbf C}f_{\mu|\mathbf c} \right)\\
&= \arg\max\limits_\mu\;\sum\limits_{\mathbf c \in \mathbf C}\log\;f_{\mu|\mathbf c}\text,\end{align}(the logarithm might be very handy, considering the exponential nature your signal of interest seems to have). From there, it's an optimization problem.


Trying to understand WAVELET TRANSFORM Frequncy-Time diagram I am reading

The Wavelet Tutorial

Part III
MULTIRESOLUTION ANALYSIS
&
THE CONTINUOUS WAVELET TRANSFORM

by
Robi Polikar

The author explains about the following fig and says;

  Note that boxes have a certain non-zero area, which implies that the
  value of a particular point in the time-frequency plane cannot be
  known. All the points in the time-frequency plane that falls into a
  box is represented by one value of the WT.

I am unable to interpret this fig and also that why value of a particular point in the time-frequency plane cannot be known.   I can always intersect two lines one passing through a particular value of t and another through f. The intersecting point will give me exact location of frequency at particular time (and V.V). I am not able to understand what thsi diagram is all about and what is its significance?

I believe by the WT, you are talking about the discrete wavelet transform, DWT.

This can be thought of as a subsampling of the continuous wavelet transform, CWT. In the case of the DWT, we pick frequencies of the form $2^{j-1}$ for ($j=1,2,\dots$) and then pick times seperated by multiples of $2^j$.

You can see this in the diagram: As frequency increases, boxes double in height and half in width.

I believe what the author means by the statement

Note that boxes have a certain non-zero area, which implies that the value of a particular
point in the time-frequency plane cannot be known.

is that the DWT, is an 'averaging' of the CWT within the boxes i.e. the boxes represent only one value and that is the average of the the values within it from the CWT. Hence as the boxes have non-zero area, it is impossible to know any of the exact values making up this average, so it is impossible to know the exact value of any point in the time-frequency plane.


How to derive the bin width / maximum frequency of DFT/STFT The bin width is given by the $F_s / N$ and the maximum frequency is given by $F_s / 2$ where $F_s$ is the sample rate and $N$ is the number of samples of the DFT.

Going between these is trivial, but I'm having trouble seeing why these formulas are correct.

I see that these formulas imply that for this image:

  halving the window duration would remove every second frequency bin, but keeping the same window duration and removing every second sample (halving the sample rate) would correspond to removing the right half of the bins.

Edit

Thank's to Dan's answer I now understand that the DTF is $N$ bins periodic ($N/2$ bins periodic for a real signal). But I want to see why either the bin width is $1/T$ where $T$ is the duration the samples are taken over, or why the DTF is $F_s$ Hz periodic. (I think bin width as $F_s / N$ is a red herring)

Aside
Sam informed me that my question assumes that the DFT size is equal to the window length and this is just a special case.

The maximum frequency for a real waveform is $F_s/2$, while the maximum frequency for a complex waveform is $F_s$.  The DFT is a transformation of N samples in time to N samples in frequency, with bins numbers $0$ to $N-1$, with bin $0$ representing "DC", and bin $N-1$ representing 1 bin less than the sampling rate.

Complex signals are resented as phasors rotating in time in either a counterclockwise ($e^{j\omega t}$) or clockwise ($e^{-j\omega t}$) direction representing positive and negative frequencies specifically. Know that each bin in the DFT represents a single complex frequency ($e^{j\omega t}$), and given Euler's Identity you need two complex frequencies to represent a sinusoidal tone:

$$2\cos(\omega t) = e^{j\omega t} + e^{-j\omega t}$$

Hence the reason why you only need $F_s/2$ for real waveforms since the second tone is redundant (always the complex conjugate for real signals).

Also importantly, the DFT is cyclical as it is a sampled system, no different than how aliasing occurs when we use an A/D converter.

The simplest way to see all of this is to consider a bicycle wheel representing a complex tone, and a strobe light as our sampling clock.

The wheel can spin counter-clockwise or clock-wise, and we will refer to counter-clockwise as a negative frequency, and clockwise as a positive frequency.

If the wheel is spinning 2 cycles per second clockwise, and our strobelight is going at 10 cycles per second, we will witness the 2 cycles per second spin rate. However if our wheel is spinning 12 cycles per second clockwise, it will appear just as if it is spinning 2 cycles per second.  If our wheel was spinning 10 cylces per second, it would appear as it is not spinning at all (DC!).

This is just like the DFT, and why bin N-1 is the highest frequency. The next bin N (if we had more bins) would be the sampling frequency itself and would be no different than bin 0. This is the cyclical nature of the DFT, and how for complex signals we can map them either from $-F_s/2$ to $+F_s/2$ or as is done in the DFT from $0$ to $F_s$. (since it is cyclical either can be a unique span).

So we have N bins with each frequency as $0/F_s, 1/F_s, 2/Fs ... (N-1)/F_s$, we see that the bin width is $F_s/N$.

Ultimately and more universally a time sample of duration $T$ has a frequency resolution of $1/T$ (Time Frequency Duality). So if our sampling rate is $F_s$ samples / sec and we have $M$ samples, then the duration $T$ will be $N/F_s$ and then our frequency resolution is the inverse of this which is $F_s/M$. If we zero pad that out to $N$ total samples, the frequency spacing of each bin will decrease to $F_s/N$ but we have not done anything to improve the frequency resolution, or our ability to discern two closely spaced frequencies- it will still be $F_s/M$! So here we note the interesting point that our frequency resolution is really not to do with the sampling rate but to do with the total duration in time of our actual (non-zero padded) data set.  Here $T$ is the length of time for the actual data within the window. With zero padding we increase the number of DFT bins which therefore increases the number of frequency samples but it does not increase the frequency resolution!  See this post for more details on that with a specific example:

What happens when N increases in N-point DFT

This assumes a rectangular window, which has the narrowest frequency resolution as given by the relationship above (and specifically this is the equivalent noise bandwidth). Any other window used will decrease the frequency resolution further, as detailed in many other posts on here with regards to windowing specifically.

Here is an example demonstrating real and complex frequencies and their DFT:

N = 10;        # Total number of samples
n = 0:9;       # sample index
k = 3;         # frequency index (Fs = N)
y = cos(2*pi*k*n/N);
stem(n,fft(y))

The plot above is for the DFT of the function $y = \cos(3\omega_n)$ where $\omega_n$ is $2\pi/N$. Note that this is the same as $\frac{1}{2}e^{j3\omega_n} +\frac{1}{2}e^{-j3\omega_n}$. Because of the cyclical nature of the DFT described above, $e^{-j3\omega_n} = e^{j(10-3)\omega_n} = e^{j7\omega_n}$, so we see the two DFT bins representing this at bin 3 and bin 7.

Think of the bicylce wheel and the strobe light, if it was rotating 7 cycles per second in the positive direction with the strobe light at 10 cycles per second, it would visually appear to be rotating 3 cycles per second in the negative direction! The upper bins of the DFT represent the negative frequencies. For real signals the negative frequencies are the complex conjugate (equal magnitude, opposite phase) of the positive frequencies, hence we only need to show the frequencies from $0$ to $F_s/2$ in the case of real signals, as we already know what the negative frequencies are.

Below shows the same result for a complex tone at at bin 3 only.

y2 = exp(j*2*pi*k*n/N);
stem(n,fft(y2))

As for the statement by the OP:

  halving the window duration would remove every second frequency bin,
  but keeping the same window duration and removing every second sample
  (halving the sample rate) would correspond to removing the right half
  of the bins.

Halving the DFT duration (not the window duration but the window duration will be modified if the DFT duration is less than double the window) would remove every second frequency bin and not impact the remaining bins only under the condition the condition that the contents remaining in the window are either cyclically invariant to the original data set (cascading the remainng window with itself results in the same original data) or there is nothing in the portion of the data that was removed (removing zeros only).  Please this post which details this further, knowing what is shown in frequency at that post can equally be described as the time domain.

Effect of changing sample rate, window duration and zero padding on DTFT and DFT


Can two nonzero signals $x[n]$ and $y[n]$ give a zero convolution Suppose $x[n]$ and $y[n]$ are two nonzero signals(i.e., $x[n] \neq 0$ for at least one value of n and similarly for $y[n]$).Can the convolution between  $x[n]$ and $y[n]$ result in an identically zero signal? In other words, is it possible that $\displaystyle\sum_{k = -\infty}^{k = +\infty}x[k]y[n-k] = 0$ for all n.

Yes, for example let

$$x[k]=1$$

for all $k$ and

$$y[k] = \begin{cases}1 & k=0\\-1 & k=1\\0 & otherwise \end{cases}$$

It is easy to see that in case of a convolution, the result will be zero for all values of $n$.


Calculating total power for the signal In my previous question I asked about calculating a single number as "Spectral Density" feature from signal data. We concluded that it is really a total power and article (link) authors wrote the wrong feature name.
I have a signal with 1 measurement per minute for a few days. It looks like this (in time domain obviously):
,timestamp,activity
0,2003-05-07 12:00:00,0
1,2003-05-07 12:01:00,143
2,2003-05-07 12:02:00,0
3,2003-05-07 12:03:00,20
4,2003-05-07 12:04:00,166

I can assume the noise is white and also other assumptions simplifying the calculations if needed (authors weren't exactly advanced in signal processing, they only used very simple features).
Questions:

How can I calculate the total power? It is "the variance plus the mean squared", so I should just calculate $(\sum_{i}^N x[i]) / N$ for the signal of length $N$? Or just sum the PSD (like in this link)?

Should I perform this calculation in the time domain (just time series like above), or should I calculate FFT and PSD (Power Spectral Density) first and calculate on PSD (which is also just an array of numbers, if I understand correctly)?

I have seen quite a few formulas for the total power. For finite number of samples and discrete measurements I have seen $(\sum_{i}^N x[i])$ divided by $N$, $N+1$, $2N$ or $2N+1$. Which one should I use? Does it have something to do with bias/unbiased sampling, i.e. in statistics we add +1 to the divisor to get unbiased calculations?

EDIT: above there should be $x^2[n]$ instead of $x[n]$

How can I calculate the total power? It is "the variance plus the mean squared", so I should just calculate $(\sum_{i}^N x[i]) / N$ for the signal of length $N$?

$(\sum_{i}^N x[i]\color{red}{^2}) / N$ is the average power, not the total energy. Just the sum $\sum_{i}^N x[i]\color{red}{^2}$ is what you're looking for.

Or just sum the PSD (like in this link)?

and

Should I perform this calculation in the time domain (just time series like above), or should I calculate FFT and PSD (Power Spectral Density) first and calculate on PSD (which is also just an array of numbers, if I understand correctly)?

The sum of the power in time domain (as above) and in frequency domain are the same, so that's just another way to compute the same (probably just more labor-intense).

I have seen quite a few formulas for the total power. For finite number of samples and discrete measurements I have seen (∑Nix[i]) divided by N, N+1, 2N or 2N+1. Which one should I use? Does it have something to do with bias/unbiased sampling, i.e. in statistics we add +1 to the divisor to get unbiased calculations?

Since you're really not looking for the power but the energy, neither. You don't have to divide at all.


Gardner timing recovery and constellation diagram I would like to ask for some guidance regarding 2 issues mainly QPSK and Gardner algorithm for timing recovery.

I am having as input a set of IQ samples (10 samples per symbol) to which I need to provide timing recovery.

Having red Gardner’s landmark paper on that concept (A BPSK/QPSK Timing-Error Detector for Sampled Receivers) I understand that in a simple BPSK case where transition exists (in a pulse always) the two first terms (x[nT]-x[(n-1)T]) of the formula e={x[nT]-x[(n-1)T] }x[nT+T/2] cancel out and the error approaches zero. In cases where no transition exists then we don’t get the desirable cancelation but statistically that is going to be resolved later due to negative symbols inserted.

All that above is quite straightforward the complications start with the obscurity behind the input samples I am having and the relation they have with the constellation diagram and the initial pulses.

Above you can see the constellation diagram (of my samples) plotted. The axes are scaled this way because the numbers are represented in 8 bit integers.
Now if I choose a smaller chunk of my sample set let say 300 and scale them properly (divide by 256) I get that :(only for clarification purposes to show the samples)

Is apparent by the concept of timing recovery itself that needs to be kept only,a number of samples that fall very close to the 4 constellation points and discard all the rest.

Here I am having two main doubts the first is fundamental related to QPSK.
How do I connect the constellation diagram with the actual waveform?

Gardner was explaining his concept utilizing 2 samples only.In my case that 10 samples per symbol are utilised means that for each pulse am sampling 10 times (sorry for the oversimplification) and the I  part is the magnitude and Q the angle of each sample?

Is it like splitting the initial wave or pulse in 10 samples (sorry again for that simplifications I know they sound trivial to many but are necessary to clear my doubts) instead of 2 that Gardner was explaining?

If yes then how do I employ the formula,I  mean, how do I set the middle  (x[n*T+T/2])  term and how the other two x[nT]-x[(n-1)T]?

However, In QPSK things are more complicated due to the fact that both I and Q samples exist.

That introduces a new term in the Gardner’s equation e={x[nT]-x[(n-1)T] }x[nT+T/2]+{y[nT]-y[(n-1)T] }y[nT+T/2]  .Where, the second term refers to the Q samples,and introduces a greater difficulty on realizing where and how to used the formula.

Any ideas on how to employ that for QPSK as well?

Thanks a lot

George

Your timing error detector should be run once for each input symbol. As you noted, your signal is far more oversampled than you need to be, so you won't end up using most of them for the purposes of timing error detection. As the paper suggests, the detector uses three samples, each separated by half the symbol period $T$.

So, for the first timing error calculation, you would use $x[0], x[5], x[10]$, for the next you would use $x[10], x[15], x[20]$, and so on. Over time you will develop an estimate of the timing error, which you can use to select the correct sampling phase for the soft symbol decision. This will reduce your signal down to one sample per symbol, giving you a (hopefully) QPSK-looking constellation.

While you aren't using most of the samples in your signal for timing error detection, they could still be useful to you during the selection of the appropriate sampling phase. This will in general require sub-sample precision, so the increased sample rate can help you in interpolating to find the signal's value at the exact position that you want.


Is there an algorithm to compare several pictures? I'm looking for an algorithm to store small dimensions of some pictures and later search on the stored data for matches similar to a new picture. Maybe something with Self Organizing Map in machine learning ?

It looks like you are looking for algorithms for Content-Based Image Retrieval. As long as you don't expect any miracles, you should be quite happy with a colour histogram approach. Have a look at MPEG-7 colour descriptors.


Why is FFT-based convolution efficient only for signals of large size? According to the documentation of scipy.signal.fftconvolve

This is generally much faster than convolve for large arrays (n >
~500), but can be slower when only a few output values are needed, and
can only output float arrays (int or object array inputs will be cast
to float).

Why does fft-based convolution is efficient only for large signals and how large they have to be (why the choice n > 500 is used)?

The direct convolution approach has a complexity proportional to $n^2$
The FFT based approach has complexity $n \log(n)$.
Since there are unknown pre-factors and the complexity only holds asymptotically, the only thing one can infer from the complexities is that only for large enough $n$ the FFT approach becomes more efficient.
The $n=500$ threshold is thus probably determined empirically.


Transfer function from experimental time domain data I have two sets of time domain data.  One is the input to a system and the other is its corresponding output, both measured at the same sampling frequency.  How to calculate the system's transfer function from this experimental data?

assuming your sets of time-domain data begin with $n=0$, Z transform both the input data, $x[n]$ and the output data, $y[n]$.

$$ X(z) \triangleq \mathcal{Z}\{x[n]\} = \sum\limits_{n=0}^{\infty} x[n] \ z^{-n} $$
$$ Y(z) \triangleq \mathcal{Z}\{y[n]\} = \sum\limits_{n=0}^{\infty} y[n] \ z^{-n} $$

divide the Z transform of the output by the Z transform of the input.  the result of that division is your transfer function.

$$ H(z) \triangleq \frac{Y(z)}{X(z)} = \frac{\sum\limits_{n=0}^{\infty} y[n] \ z^{-n} }{\sum\limits_{n=0}^{\infty} x[n] \ z^{-n} } $$

i cannot guarantee that, if you were to factor numerator and denominator, there will not be identical poles with zeros that could cancel and make the transfer function simpler, but it's a transfer function.


Fir filtering operation? Also convolution? While reading signal processing first

I came across a concept which indicated that fir filtering operation is also convolution
As shown in attached photo

In figure caption line,convolution is written in bracket after fir filtering

If yes?then what about IIR filtering operation? Is that also convolution?

Yes in general both IIR and FIR filters can be derived from the convolution summation, but since and IIR by definition has an impulse response that goes to infinity, to implement this directly would require an infinitely long convolution. The FIR can be implemented from the convolution expression directly, but a recursive subset from difference equations is used to implement IIR filters.  So convolution still takes place in both cases (the output is the convolution of the input with the impulse response), except in the case of the IIR filter the generalized convolution expression cannot be used directly for implementation.

Given the generalized convolution summation as:

$$y[n] = x[n] \star h[n] = \sum_{k=-\infty}^\infty h[k]x[n-k]$$

$h[n]$ is the impulse response of the filter and $x[n]$ is our input signal.

Here is a convolution with a Finite Impulse Response (FIR), which is causal (starts at $k = 0$, and has a finite impulse response (summation is over $M$ samples) As written this is directly realizable as an FIR filter:

$$y[n] = \sum_{k=0}^{M-1} h[k]x[n-k]$$

In comparison here is a convolution with an Inifinite Impulse Response (IIR), which is also causal but the impulse response extends to infinity:

$$y[n] = \sum_{k=0}^{\infty} h[k]x[n-k]$$

This is impossible to create a direct realization, but a subclass of IIR systems can be implemented using recursive difference equations:

$$y[n] = x[n] -  \sum_{k=0}^{N} h[k]x[n-k]$$


Downsampling in MATLAB - why scaling filter coefficients? Trying to undertsand how MATLAB does resampling, looking at toolbox/signal/resample.m.

This is the part where (downsampling) filter coefficients are calculated:

  fc = 1/2/pqmax;
  L = 2*N*pqmax + 1;
  h = firls( L-1, [0 2*fc 2*fc 1], [1 1 0 0]).*kaiser(L,bta)' ;
  h = p*h/sum(h);

(In this case pqmax and p both represent downsampling ratio.)

Not sure what is the purpose of last line where all filter coefficients are scaled by downsampling ratio over coeffs sum, p/sum(h)? Does anyone know the theory behind?

MATLAB firls
MATLAB kaiser
firls theory: Least Squared Error Design of FIR Filters

This is done to achieve unity gain, i.e. to have the same energy in the downsampled and original signal.

The sum of filter coefficients sum(h) is the filter gain, by dividing the filter by it, it is set to 1 (unity).
It is multiplied by p to compensate for the missing samples, again to achieve the same engery as in the original signal.


FIR Filters and moving average filter Fir filter confusion

I am trying to get some knowledge about FIR filters. According to above link FIR filter is just the convolution of a sampled input signal to impulse response of a filter. I have done some research and according to some texts FIR is the same as calculating the moving average of the input samples and there is no talk about impulse response. It seems if we take the moving average we can create a filter.

I would appreciate it if someone could help me to understand FIR filter . According to tutorials here we make FIR filter by just delaying the samples .https://greatscottgadgets.com/sdr/10/

A FIR filter is (similar to, or) just a weighted moving average filter, where (unlike a typical equally weighted moving average filter) the weights of each delay tap are not constrained to be identical or even of the same sign.  By changing various values in the array of weights (the impulse response, or time shifted and sampled version of the same), the frequency response of a FIR filter can be completely changed.


Phase function of filter I have a filter with the transfer function

$$H(z) = 1 - 2z^{-2} + z^{-4}.$$ The task is to find the phase function $\theta (\omega).$

My attempt is to start by expressing the frequency response
\begin{align}
H(\omega) &=H(z)\big\vert_{z=e^{j\omega}}\\
&= 1 - 2e^{-2j\omega} + e^{-4j\omega}\\
&=  e^{-2j\omega}( e^{2j\omega} - 2 + e^{-2j\omega})\\
&=  e^{-2j\omega}(2\cos(2\omega)-2)
\end{align}

I think there is a relation of the sort $$H(\omega) = e^{j\theta (\omega)}|H(\omega)|$$

which in my case would give me $$\theta (\omega) = -2\omega.$$

The correct answer is $$\theta (\omega) = -2\omega + \pi.$$ How do I get the extra term?

HINT:

Answer these two basic questions, and you should be able to solve your problem:

What's the sign of the term $2\cos(2\omega)-2$? Can it be a magnitude?
What does a term of $\pm \pi$ in the phase do to the sign of an expression?


Fourier Transform of triangle function $x(t)=\Delta\left(\frac{t-1}{2}\right)$ Can you please tell me if my working is right for the Fourier Transform of this function: $$x(t)=\Delta\left(\frac{t-1}{2}\right)$$

My workings are:

I have used the fourier transform standard results. Please could you tell me if I which I have to apply:

First the $F(t-1)$, then $F(t/2)$
First the $F(t/2)$, then $F(t-1)$
First the $F(t/2)$, then $F(-1/2)$

Please tell me where I have gone wrong.

My final answer is:

$$x(t)=\mathrm{sinc}^2(\omega/2)e^{-j\pi f}$$

I will not go into the details of the specific function that you need to transform. This is for two reasons: first, with a bit of effort you will be able to do it yourself (after having read this answer), and second, because you didn't define your functions. Note that it is not immediately clear how $\Delta(t)$ is defined; also note that there are two frequently used definitions of $\text{sinc}(x)$ (either with or without an implicit factor of $\pi$).

I will discuss the way to derive the Fourier transform of a shifted and scaled version of a function $f(t)$ with Fourier transform $F(\omega)$. For this you need to know the two following Fourier transform properties (derive them yourself and you'll never forget them!):

$$f(at)\Longleftrightarrow \frac{1}{|a|}F\left(\frac{\omega}{a}\right)\\
f(t-a)\Longleftrightarrow F(\omega)e^{-ja\omega}\tag{1}$$

The Fourier transform of $f((t-a)/b)$ can be derived from $(1)$ in two different ways:

Find the transform of $f(t-\frac{a}{b})$, then substitute $t\rightarrow t/b$
Find the transform of $f(\frac{t}{b})$, then substitute $t\rightarrow t-a$

The first way is then
$$f\left(t-\frac{a}{b}\right)\Longleftrightarrow F(\omega)e^{-ja\omega/b}\\
f\left(\frac{t-a}{b}\right)\Longleftrightarrow |b|F(b\omega)e^{-ja\omega}$$

The second way of course gives the same result:
$$f\left(\frac{t}{b}\right)\Longleftrightarrow |b|F(b\omega)\\
f\left(\frac{t-a}{b}\right)\Longleftrightarrow |b|F(b\omega)e^{-ja\omega}$$


Can IR emitter signal be distorted by curved glass housing around receiver? I have built a garage door opening circuit that includes an IR sensor and an IR emitter (a remote control). When the IR sensor is in the open air it works fine.

If I put the sensor inside a glass jam jar, the signals received are different but sometimes consistent, ie. each button on the remote returns a consistent (wrong) signal from the receiver.

My understanding is that the signal is a combination of high/low pulses with different duration.

Could the jam jar's reflection and refraction be impacting the length of the pulses? If not, how can this effect be explained?

Short answer: No.

Long answer: you can of course create shadowing that way, and that would disturb operation. And of course, a glass wall will refract infrared just as it refracts any other light. But glass is usually IR-transparent, just as it is transparent visible light, and an IR remote isn't exactly a source of a focused beam, so chances are that, no, there's no difference in operation.

What will happen is that the IR will be a tiny bit weaker, due to absorption in the glass. But that won't be really much, I guess.

  Could the jam jar's reflection and refraction be impacting the length of the pulses?

How so? It's the same signal entering and exiting the glass. Admittedly, at a different position, maybe a bit attenuated, but still, if you wanted to change the duration of the pulses, then you'd need a frequency converter (which is a nonlinear device) or a time machine. You have neither.

What could be happening is that due to reflections in the glass, at the receiver, you'd get one received signal that took "the shortest" path, overlaying with a "delayed" path. For that to have any negative effect, the path length difference would need to be in the order of the duration of your pulses.

Infrared remotes work with longer pulses modulated to something like 38 kHz, if I remember correctly. Ignoring that wavelengths are different in glass (factor of maybe 1.5, and we're just here to get a feeling for the orders of magnitudes involved): a signal travelling a speed of light travels, in one period of a 38 kHz signal,
$$\frac{3\cdot 10^8\frac{\text m}{\text s}}{3.8\cdot 10^4\,\text{Hz}} \approx 0.7 \cdot 10^4\,\text m = 7\,\text{km}$$

Since any distance your "detouring" beam might take will be way, way smaller than that, we can rule that out.

I think for the typical roughness of a marmelade jar glass, we can also rule out thin-film effects on the surface.


How do I analyze a music file ( only contains 5 seconds of note A) using the FFT plot to determine the note being played? To perform this task why can't I just import the file, and directly perform the fft and use the stem plot to determine the note? From the index k of the peak I can do w = 2pi*k/N where N is total number of points.

Fs, x = wavfile.read('piano.wav')
X = np.fft.fft(x)
plt.plot(abs(X))
plt.show()

However I'm not sure about the correct way to map the fourier coefficients from the plot to the frequency to determine the note.

Please advise.
Thanks.

An FFT magnitude peak frequency (derived from a single FFT result index) is not reliably the same thing as the music pitch frequency or the musical note.  Note pitch is a psychoacoustic phenomena, and more reliably related to periodicity than the (sinusoidal) spectral peak frequency (e.g. a periodic/repeating waveform shape does not need to look like a sine wave).  See FFT Pitch Detection methods: Autocorrelation or other? for a list of some better methods for estimating note pitch frequency or the musical note.


Maximum dB SPL from a digital microphone Basic sanity-check question: if I have a digital microphone with a sensitivity of -26 dBFS (say, this one), does this mean that the max dB SPL it can handle is 120 dB SPL?

My thinking is that if the sensitivity is measured at 94 dB SPL and the max output is presumably 0 dBFS, then that 26 dBFS sensitivity would bring us up to 120 dB SPL. Is that right?

Mostly yes. See http://www.analog.com/en/analog-dialogue/articles/understanding-microphone-sensitivity.html

You need to be aware of the difference between peak and RMS level. The max SPL for a sine wave would indeed by 120 dB SPL. However the max SPL for something like a trumpet would be significantly lower, at least if you define SPL as "what I would see on an SPL meter.

The difference is the crest factor of the wave form: for a sine wave with an RMS of 120 dB, the peak level is 123 dB. For a trumpet with 120 dB RMS, the peak level will be around 140 dB or so. So the microphone will clip way earlier on a trumpet as compared to a sine wave.


How to detect quadrilaterals in a point cloud? How can I detect a quadrilateral in a point cloud that looks almost like a rectangle? How do I set the criteria (i.e orientation, size etc) to get the shape for example in the image below?

There is one possible solution which is roughly $O(n\cdot (n-1)\cdot log(n))$ where $n$ is the number of points.

Build 2-dimensional k-d tree from all the points
For each pair of unmarked (undiscovered) points

2.1 Determine two areas (based on rectangle criteria) where the other two points can possibly lay

2.2 Find nearest point to center of each of the areas (using the k-d tree)

2.3 If each nearest point is unmarked and lays within its respective area, store new rectangle, mark all four points as discovered

The algorithm can be optimized by keeping undiscovered points in a set, picking new ones from the, progressively reducing its size.

Another speedup would be removing discovered points from the KD-tree so the algorithm does not have to check them in step 2.3.


From fft result one can affirm that the signal is affected by white noise? I calculated the FFT of the Az acceleration sampled with MPU6050 connected to Arduino.

Wikipedia definition: "White noise" is  a random signal having equal intensity at different frequencies, giving it a constant power spectral density

Can I say from my FFT that the reason of the variance of the signal, in my test condition, is caused by white noise?

Additional question on FFT

If my sampling frequency is 56Hz are my bin_frequency 28?
My frequency bin is influenced by the length of the input vector or only by my sampling frequency?

Test condition:

I recorded of the z value in static condition on a table, without any perturbation
The signal was  nonuniformly sampled so I resampled it

Input

Sampling frequency 55.556

Az 718X1 Vector

Std_Input_Vector_Resampled 0,0385

FFT Output

FFT 718X1 Vector

Without Zoom

ZoomX3

ZoomX4

This is my first time I approach the FFT(I've been studying this argument from one month so I am very noob), and I have a civil engineer background, so is my first time with signal analysis, never studied Signal Theory before.

My interest for the FFT is to define the best filter to identify the pothole and road anomalies

Thanks for your time, thanks for your patience.

These are just some hints on how to proceed which are too long for a comment. We cannot give a definite answer from just looking at the plots.

First, you should plot the magnitude of the FFT in dB, i.e. plot(20*log10(fft(signal))). Then, you would not need to plot different zoomed versions.

To refer to question 1) Is this the whole signal or just a part of it? Try analyzing a longer measurement cycle to become more certain. Certainly, there is some addition of white noise to the signal. But, additionally, it seems to me there is an additional sine-wave of period length 4 contained (referring to your first plot).

To do some more analysis:

Assuming you measure a constant acceleration, remove the mean from the signal. Then, the remaining signal just contains the measurement error.
Have a look at the (circular) autocorrelation of this error. In case of white noise, it should be ideally a single peak at delay 0. If the autocorrelation exhibits extra (significant) peaks, there probably is something extra on top of the white noise. (Looking at your FFT-plots, I actually expect some extra peaks).

To answer questions 2), 3) (I dont really understand your question 2)):

If your sampling frequency is 56Hz, your FFT output will range from 0 to 56Hz, where the distance between the bins is 56Hz/FFT-length. Furthermore, as your input signal is real, the frequencies above 56/2Hz will be a mirror of the first half of the frequency domain output, they do not contain additional information. As a general rule for Fourier Analysis: The frequency resolution $\Delta_f$ after FFT is only determined by the time duration $T$ of the original signal, given by $\Delta_f=1/T$. You can have a look at this article I wrote on my website about the FFT and frequency domain resolution as an introduction to the topic.


High-pass filter with quantized (fixed-point) coefficients I have a high pass IIR filter with numerator
b = [0.997688673703782,-1.995377347407564,0.997688673703782];

and denominator
a = [1,-1.995372005171182,0.995382689643947];

the frequency response seems ok with fvtool(b, a). However, the response changes a lot after the coefficients are converted to fixed point numbers by
bb = round(b * 32768) / 32768;
aa = round(a * 32768) / 32768;

The frequency responses are as follows

How do I design a fixed-point digital filter with specific low frequency response? For example, a high pass filter with a desired cut off normalized frequency. Thx!

One problem with your design is that the quantized numerator coefficients don't add up to zero anymore, so you lose the desired notch at DC. That can be done differently, because even with quantized coefficients you can have a (double) zero at DC. Just make sure that $b[0]=b[2]$ and $b[1]=-2b[0]$.
The real problem is the denominator polynomial. It is well known that the available pole locations resulting from quantizing the polynomial coefficients are very sparse at low frequencies.
Note that the maximum magnitude of the ideal coefficients is close to $2$, so in order to quantize the coefficients to $16$ bits, you need to scale by $2^{14}$ instead of by $2^{15}$. This will result in the correct value range of the quantized integer coefficients:

bb = round( b * 2^14 ) * 2^(-14);
aa = round( a * 2^14 ) * 2^(-14);

The resulting quantized numerator coefficients guarantee a double zero at DC. However, the quantized denominator coefficients result in two real-valued poles, one of which is on the unit circle. I.e., there is a pole-zero cancellation and you end up with a first-order filter.
So one of the options you have is to use a first-order filter instead of a second-order filter. The numerator coefficients are trivial in this case:

b1 = [1,-1];

The denominator coefficient can be chosen experimentally using some negative value close to (but smaller in magnitude than) $-1$. Of course that value must be representable using $16$ bits:

a1 = [1, round( -r * 2^(15) ) * 2^(-15)];

The plot below shows the frequency response of the ideal second-order filter (blue) and of the quantized first-order filter with $r=0.998$. You have to judge for yourself if the deviation from the ideal response is tolerable for you application or not.

Another solution is to use a different filter structure which is more robust to coefficient quantization (at least at low frequencies). Such a structure is the coupled-form structure, which is described in most textbooks on DSP (e.g., in Oppenheim's Discrete-Time Signal Processing).


Discrete Fourier transform of complex time series I am working through previous exam questions for a class I'm currently taking in signal theory, and one problem, although it looks easy, stumps me.

PROBLEM

Two real $N$-periodic time seires $f_n$ and $g_n$ are given, where $n = 0, 1, 2, ..., N-1$.

In addition, the complex time series $h_n$ is given as

$$h_n = f_n + i g_n$$

and the discrete Fourier transform

$$H_m = F_m + i G_m$$,

where $m = 0,1,. . ., N-1$.

Find the discrete Fourier transofrm $F_m$ and $G_m$ to the time series $f_n$ and $g_n$ expressed through $H_m$

ATTEMPTED SOLUTION

As I said, I'm really stumped here.  Other than writing:

$$F_m = H_m - i G_m$$

and

$$G_m = i (F_m - H_m)$$

I don't see what else to do.  Is this all there is to it?  It seems way too easy, so I'm pretty sure I'm not doing this correctly.

If anyone can give me any help, I would greatly appreciate it!

Proposed answer to my own question:

We have that:

$$F_m = \frac{1}{2}(H_m + H^*_m)$$

and

$$G_m = \frac{1}{2i}(H_m - H^*_m)$$

where $H^*_m$ represents the complex conjugate of $H_m$.


How to add noise to time series data in Python? I have time series data with $x$, $y$, $z$ sampled at $50\textrm{ Hz}$ from a sensor and I would like to add a uniform and Laplace noise to it. How would I achieve this in Python?

Some comments:

I'm not sure if this is the right forum for python questions.
You left some relevant information, like the type of variable your data is stored.
I assume that "la place" is laplace distribution.

Please see if this helps, and make the necessary adjustments:

import numpy as np
def fn_addnoise(data):
    i = len(data)
    # create 1D numpy data:
    npdata = np.asarray(data).reshape((i))
    # add uniform noise:
    u = npdata + np.random.uniform(size=npdata.shape)
    # add laplace noise:
    p = npdata + np.random.laplace(loc=0.0, scale=1.0, size=npdata.shape)

    print npdata
    print u
    print p
    return u,p
# 'data' might be a tuple...
data = (1,2,3,4,5,6)
# or a list...
data = [1,2,3,4,5,6,7]
fn_addnoise(data)


Audio texture analysis for visualization I'm making a complex music visualizer as a personal project. I've already determined that I will be using a series of modified Goertzel algorithms for frequency detection to place the frequencies on a spiraling piano scale without losing all detail in the bass or calculating unnecessary harmonics.

I'm now looking for analysis for the post-processing effects. I'm looking to estimate the warmth/brilliance (or liveness), and 'grittiness' of the sound as I go. A good example of what I mean by grit is "Dark, Darker, Yet Darker."

How should I approach this problem? I don't have a large, supervised dataset and I only have access to a laptop at the moment, so training a neural net is pretty much out of the question for me. For the record, I'm using Processing (java-based) as my coding environment.

Listening to that 'gritty' track you posted, two qualities stood out to me: firstly the bitcrushing / 8-bit effect, and secondly the overdriven quality of the guitar. Overall the effect sounded like there was all this 'noise' in the background, and the pitches of the instruments were competing with this.

Some thoughts based on that observation:

If you're already doing frequency detection, maybe you could you measure the proportion of energy at/near these frequencies, as a proportion of the total energy. This is almost like an estimate of SNR. So based on that, you could say that a higher 'SNR' is gritty sounding.
Alternatively, you could do some sort of frequency domain correlation with an ideal noise spectrum of some sort, e.g. white or pink noise.
Based on this page, both bitcrushing and overdrive create square waves in the waveform. So you could try and measure how square the waves are, perhaps by correlating the waveform with a superposition of square waves at your detected frequencies.


Automatically detect an image collage I'm trying to automatically detect if an image is a collage vs a single photograph. I'm not too concerned with difficult cases of artistic collages. What I'm trying to solve is rectangular collages like below. I've tried edge detection (canny) + vertical and horizontal sobel filtering + line detection (Hough transform) to try to identify perpendicular lines but am getting too many false positives. I'm not very good at image processing so any input would be welcome. Thx!

Here is a basic answer with some rapid experiments to show you what is possible to do with Morpholocical operations. Note that it is just a proof of concept, not a viable solution and you should improve it to better fit your problem.

First filter your image to detect edges. I used Canny but you may use a more appropriate one.

Then filter the BW result to keep only vertical and horizontal (hypothetical) lines by eroding the result with line structural elements. Here is the result for the erosion with a horizontal line of 5px:

Lastly sum each eroded images on its corresponding dimension. Here is the result of the sum of the previous image on the 2nd dimension:

From here you have to chose a threshold in order to infer the presence (or absence) of a collage. In the previous graph the sharp peaks indicate the position of horizontal boudaries of frames. A threshold value of 100 may be appropriate as it indicated that more that $1/3$ of the horizontal line is populated with edge pixels.


Running Integral of sine and cosine functions In typical signal processing course we were taught that the integral of signal $x(t)$ is given by
$$y(t) = \int_{-\infty}^{t}x(\tau) d\tau$$
How can we use this definition to evaluate the integrals of eternal sinusoids, $x(t) = \sin (t)$ or $x(t) =\cos(t)$. I mean how do you evaluate / prove that $$\int_{-\infty}^{t}\cos(\tau)d\tau = \sin(t)$$

Thanks in advance.

Edit: I raised this question because of its (ab)use in frequency modulation(FM) for single tone signal. To give details, any FM wave of single tone signal $m(t) = \cos(\omega_0 t)$ is given by
\begin{align}
\Phi_{FM}(t) & = A\cos\left(\omega_c t + K_f\int_{-\infty}^{t}\cos(\omega_0 \tau)d\tau\right)
\end{align}
It can be written as
\begin{align}
\hat{\Phi}_{FM}(t) & = A\: Re\left\{\exp\left(j\omega_c t + jK_f\int_{-\infty}^{t}\cos(\omega_0 \tau)d\tau\right)\right\}
\end{align}
or
\begin{align}
\hat{\Phi}_{FM}(t) & = A\:Re\left\{\exp\left(j\omega_c t\right) \exp \left(j K_f\int_{-\infty}^{t}\cos(\omega_0 \tau)d\tau\right)\right\}
\end{align}
In communication systems textbooks, the argument of exponential term involving cosine term is magically written as $\dfrac{K_f}{\omega_0} \sin (\omega_0 t)$ which really fascinates me.

Later they claim this is periodic signal and hence can be represented as Fourier series expansion involving Bessel functions etc.,

As I mentioned it is really magic and am interested in knowing this magic trick. Hence this post. Sorry for being so long....

Thanks to all of you.

I would not teach integrals this way.

Issues might arise form the definition of terms. In standard calculus, a (note: "a", not "the) primitive function (also called antiderivative or indefinite integral) of a continuous function $f$ is a kind of converse of the concept of derivation. SO:

  If $F$ is "a" differentiable function  whose derivative is equal to
  $f$, then $F$ is a primitive (function) of $f$. There exists many such
  functions, up to a constant $c$, because: $(F(t)+c)' = (F(t))' =f(t)$.

Here, you can check that $t\mapsto \sin(t)+ c$ is differentiable, with derivative $t\mapsto \cos(t)$.

Now (if I remember well), suppose that $f(x)$ is a continuous function, AND that it is  "integrable" over $\mathbb{R}$, then $F(t) = \int_{-\infty}^t f(\tau)d\tau\;$ exists, and is a differentiable function. THEN, it is also a primitive of $f(t)$ (among others).

But here, it is not possible, namely because $t\mapsto \cos(t)$ is not "integrable" over $\mathbb{R}$, and thus the integral expression makes no sense.

I would say that one should not define "integrals" the way you wrote it, if one wants mathematical correctness. However, for causal signals, in a lousy engineer perspective, it may be intuitive. But the cosine is not causal.

[Note: there are more general definitions of primitives]


Subtract a digital filter from another filter How to subtract a digital filter from another one, if their lengths are different. How to make the length equal.

If (and only if) the filters (either FIR or IIR) are linear, then you can subtract their coefficients term-wise, supposing they are correctly aligned (below, on the $0$ index), treating them as infinite series with trailing zeroes when the coefficients are not defined. For instance with:

$$h = [h_0, h_1, h_2]$$
and
$$g = [g_{-1},g_0, g_1, g_2,g_3]$$

since the convolution with $h$ yields the same result as the convolution with a zero-padded $h$ (and $g$)

$$h' = [\ldots,0,0,h_0, h_1, h_2,0,0,\ldots]$$
the difference filter $d$ will have (symbolically $h'-g'$):

$$d = [\ldots,0,-g_{-1},h_0-g_{0}, h_1-g_{1}, h_2-g_{2},-g_{3},0,\ldots]$$

For FIR filters (or truncated IIR), you can limit to a finite sequence bounded by the smallest and the biggest index of both filters.


Get sound from matlab and convert with FFT My question is simple but I am new to MATLAB's sound functions so I am unsure what to do. I would like to read the audio (from the microphone) then apply the fft(X,n) function to it where the n should be 256. I was applied the fft function to this code like fft(myRecording,256) but I kept getting zeros.

I use this function:

function audio_recorder = MicFFT(Fs, update_rate)

recObj = audiorecorder(Fs, 16, 2);
get(recObj);

% Record your voice for time specified by update_rate.

while 1
    recordblocking(recObj, update_rate);
    myRecording = getaudiodata(recObj);
    subplot (2, 1, 1);
    plot (myRecording);
    xlabel('Time')
    ylabel('Amplitude')
    title('Time domain representation of y(t)')

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

    subplot (2, 1, 2);
    plot(f,2*abs(Y(1:NFFT/2+1)))
    drawnow;
    axis([0 Fs/2 0 0.01])
    title('Single-Sided Amplitude Spectrum of y(t)')
    xlabel('Frequency (Hz)')
    ylabel('|Y(f)|')

end

end

Fs is the sampling frequency whilst update_rate is the length of each recording


What Is "Description Vector" in Image Processing? When we want to use classifiers like SVM, we first should extract descriptors using algorithms like SIFT.

But I have a question which might call a silly one:
Let's assume
$$
\begin {equation}
D_a=\{d_{a,1},d_{a,2},d_{a,3},...,d_{a,M}\}
\end {equation}
$$
where M is total number of descriptors and $d_{a,?}$ denotes the description vector of mth descriptors.
My question is: what the sentence "$d_{a,?}$ denotes the description vector of m-th descriptors" means? Does it mean that we have 1000 interest point in the image that each of them has ,let's say, a 128-dimension SIFT descriptors? Do we have 128000 values?

I think you have a matrix.
Each Row / Column is a descriptor vector of a point in the image.

Just like having features, let's say M features, and each point has M values corresponding to M features.

So each element in the descriptor vector is a value of specific feature for this point.

And yes, if you have M = 128 Features and 1000 points you'll get a matrix with 128,000 entries.


Why do DFT frequency buckets need to be divided by sample period? edit: This is my first post in this community. I'm sure that my downvoter had a good reason to do so, but could someone please comment and tell me how I can better format my question?

--

I have a DFT algorithm that almost makes sense. The issue is when I go to plot the frequency response. For some reason, the frequency response gets multiplied by the sampling period (see first image).

Understanding this, I can fix it by dividing my frequency response X-axis by sampling period (see second image), but I don't understand why this is an issue in the first place.

Code is posted below - can anyone explain what I am missing?
(example uses a 1 Hz sinusoid)

Here's how I think it should work:

But here's what I actually have to do to make it work:

x_t = sin( 2 * pi * t);     % Input: 1 Hz Sinusoid
s_f = zeros(N,length(x_t)); % Output: Freq spectrum

%%%% init params
    T    = 4;     % (sec) time window
    dt   = 0.1;   % (sec) sample time

%%%% init comp
        N = T / dt;          % (ct) num samples
        n = 0 : N - 1;       % (ct) bucket index vector
        t = 0 : dt : T - dt; % (sec) time vector
smpl_freq = 1 / dt;          % (Hz) sampling frequency
 freq_res = smpl_freq / N;   % (Hz) frequency resolution
        f = n * freq_res;    % (Hz) freq. bucket vector
        b = 2*pi*(f/N) .* n.';

    anlyz_fn = cos(-b) + i * sin(-b);

%%%% run DFT
for k = 1 : N
    val = 2 * sum(anlyz_fn(k,:) * x_t.');
    s_f(k) = val;
end

%%%% Plot Results
subplot(2,1,1);
plot(t, x_t);

subplot(2,1,2);
stem(f(1:N)/T , abs(s_f(1:N))/N);

I don't know why you got downvoted.  This is actually a very pertinent question having to do with normalization conventions.

The most conventional DFT formula is this:

$$ X[k] = \sum_{n=0}^{N-1} x[n] e^{-i \frac{2\pi}{N} n k } $$

The $1/N$ normalized form is

$$ X[k] = \frac{1}{N}\sum_{n=0}^{N-1} x[n] e^{-i \frac{2\pi}{N} n k } $$

And the $1/\sqrt{N}$ normalized form is

$$ X[k] = \frac{1}{\sqrt{N}}\sum_{n=0}^{N-1} x[n] e^{-i \frac{2\pi}{N} n k } $$

My very strong preference is to use the $1/N$ normalized version.  This makes the bin values work as you expected them to.  It also makes the expected bin values at whole integer frequencies bin count invariant, and leakage values comparable.

If you do not apply a $1/N$ in the forward DFT, you will need it on the inverse to get your original signal back.  If you put $1/N$ on the DFT, then the inverse needs no rescaling (and the bin values serve as coefficients for continuous Fourier series).

The $1/\sqrt{N}$ is the "proper mathematical" one as it makes the forward and inverse DFT unitary matrix multiplications.  Yeah, nobody does that.

Bottom line, I think the $1/N$ normalization should be the standard convention and I write all my blog articles using that definition.

BTW, you will get a magnitude of $1/2$ for pure whole integer tones on $1/N$ normalized DFTs.  The other half is on the other side of the circle.