texts:
2nd order Butterworth transfer function:

$$H(z) = \frac{ 0.0201z^2 + 0.0402z + 0.0201 } { z^2 - 1.5609z + 0.6413 }$$

I'm a bit uncertain about the factorization at the denominator. Here's what I have so far:

\begin{eqnarray*}
H(z) &=& \frac{ 0.0201(z^2 +2z + 1)  }{ z^2 - 1.5609z + 0.6413 }\\
 &=& \frac{ 0.0201(z + 1)(z + 1)  }{ z(z - 1.5609) + 0.6413 }
\end{eqnarray*}
Zeros : $z + 1 = 0$, zero at $z = -1$.

How can I find the poles?


texts: What should be the dimensions of a Fourier Transformed image?

texts:
I have applied Fourier Transform to the following image.

I have downloaded this image from the Internet. I re-sized (without maintaining the aspect ratio) it using MS Paint application of Win7 to make it 512x256.

Then I have used two applications to observe its Fourier Transformed appearance.

IPLab gives the following output:

ImageJ2-20160205 gives the following output:

As you can see, the first output is a 512x256 image. The second one is a 512x512 image.

Why are those outputs different?

How would they effect processing of an image?


texts: transmitting the real/imaginary parts of pass-band signal

texts:
As I know, the transmitted pass-band signal whose carrier frequency is $f_c$ can be written as follow:
$s_1(t) = \text{Re}(a(t)e^{j2πf_ct})$ , where $a(t)$ is the complex modulated base-band signal. It's clear to notice that at the receiving end,   we can recover the complex modulated symbol $a(t)$ based on the real transmitted part $s_1(t)$.
But the question which I couldn't get a sure reply for it, what's about if we transmitted the imaginary part of our signal? .. it means that  $s_2(t) = \text{Imag}(a(t)e^{j2πf_ct})$. Will we be able to recover the complex modulated symbol based on $s_2(t)$ too?  How can we prove that mathematically?
thank you in advance


texts: Minimum phase conversion with Cepstrum method, how to scale the result?

texts:
I am trying to convert a zero phase spectrum (magnidue response curve with zero phases) to a minimum phase spectrum, because I need a totally causal impulse response for FFT spectral filtering, and having an IR devoid of anticausal components allows me to implement a zero latency partitioned scheme.
By now I experimented with two similar methods based on the real cepstrum but my question is inherent to the second method.
In the first method you compute the cepstrum by taking the ifft of the log of the magnitude curve, negate the second half of the cepstrum (actually I rather use a cosine function to have a more graceful transition and less ringing artifacts), take the fft of the result and then instead of computing the exp of the resulting spectrum I simply keep the phase combining it with the original magnitude. This method works well enough despite the resulting IR doesn't always fade to zero.
In the second method you zero the second half of the cepstrum instead of negating it (here too I rather prefer a rised cosine for a more graceful transition), double the result, take the fft and then exponentiate the result.
Here is the problem. The resulting spectrum is uniformly scaled by a variable amount which seems depending on the shape of the original magnitude spectrum, so I cannot compensate by scaling by a fixed value because I don't know the relationship between the resulting scaling factor and the original magnitude spectrum, and I can't definitely have a global level which is dependent on it! I searched but I found very little information around about the min phase cepstrum method. How can I scale the final result properly? I realize this may be a very specialistic topic. Thanks


texts: How to calculate the upsampling rate (sample per symbol) for a recorded raw IQ .wav file?

texts: I am not familiar with raw IQ wave files. For tetra-1 signal, the bit rate is $36$ k bps (or $18$ k as the symbol rate). If the Fs saved in the .wav file is $48000$ Hz, what should be the SPS (sample per symbol)? Remember that if SPS is 4 for QPSK, then both I&Q components need to be upsampled by $4$. In this case, should it be calculated as $(48000/2)/18000$?


texts: understanding the result of the simulation of In-phase and quadrature part BPSK signal

texts:
I am a student and working on a hardware implementation of a transceiver and going to check an algorithm of time recovery that includes symbol synchronization and sampling clock synchronization. The proposed algorithm has a problem in realisation in hardware and I was thinking to check the algorithm and whole simulation in Matlab and probably I could to fix and improve it.
In Matlab, I would like firstly generate a test signal which will be received. Before I asked here about the time recovery algorithms and it was explained to me excellent. Thank you, guys!
My simulation in Matlab:

100 symbols with 4 samples
random noise
raised cosine filter

Code:
% create a random sequence of 1 or -1
N = 100;
bit_stream = rand(1,N)>0.5;
BPSK_symbols = 2*(bit_stream)-1;

% Filter coefficients
h=firrcos(150,.5,.2,4);

% add 4 zeros
for i = 1:length(BPSK_symbols )
   BPSK_symbols_zzeros(1) = BPSK_symbols (1);
   BPSK_symbols_zzeros(4*i)=BPSK_symbols (i);
end

% apply filter
tx_signal = conv(filter_coeff,BPSK_symbols_zzeros);
% generate noise
noise = rand;
norm=std(tx_signal)*noise;
% add noise
for j = 1:length(tx_signal)
   phi = 0.89 % pi/4
    tx_signal_w_phi(j) = tx_signal(j)*exp(sqrt(-1)*(phi_1))+0.7*norm*(randn+sqrt(-1)*randn);
end

My result : red - inphase , green -quadrature
I think the result is wrong. If the phase offset is pi/4, I have to get I-Signal the same as Q-Signal, isn't it?
As you see, the max amplitude (peaks) are different.
Could someone help me to fix this implementation, because I can't understand where I have mistaken...


texts: Why is my resampled signal cutoff when resampling with a windowed FIR filter in MATLAB?

texts:
I have an exercise to write a function which upsamples or downsamples a signal by a factor of F (integer) using a FIR filter windowed with a Kaiser window function. Exercise I had last week included us to write a function which returns an impulse response of such filter so we can use it in this exercise rather than using a built-in function from MATLAB.

The function h = firkaiser(Wc, Bt, As, Ap) creates a NF filter with central frequency Wc, Bt which is a band where the filter cuts off, and attenuations in the stopband and passband. Using a standard procedure of calculating M (length of impulse response) and windowing an ideal filter i have created this function. Here is an example of a NF FIR filter created by this function just to make sure the problem isn't in the filter itself.

h = firkaiser(0.3*pi, 0.3, 50, 0.5)

Gives these results:

Now to the problem. In the exercise description as a hint it is mentioned that the interpolaton filters used must have Wc = pi/F, Bt = 0.4*pi/F, As = 70dB and Ap = 0.85dB so we must use these values. We also must compensate the delay so that the output signal is in phase with the input signal. Here is the function i wrote:

function [y] = resamplesignal(x, sampling, F)

if (sampling(1) == 'u') %if you want to upsample
    xr = upsample(x,F);
    fir = firkaiser(pi/F, 0.4*pi/F, 70, 0.85);
    y = filter(fir, 1, xr);
    delay = mean(grpdelay(fir)); %removing the delay
    y(1:delay) = [];
    y = F*y; %upsampling lowers the amplitudes

    figure
    stem(1:F:F*length(x), x);
    hold on
    plot(y);
    pause
else
    fir = firkaiser(pi/F, 0.4*pi/F, 70, 0.85);
    xr = filter(fir, 1, x);
    delay = mean(grpdelay(fir)); %removing the delay
    xr(1:delay) = [];
    y = downsample(xr, F);

    figure
    plot(x);
    hold on
    stem(1:F:F*length(y), y);
    pause
end
close all
end

The function works as expected but the last part of the signal is lost first because of the delay and when the delay is deleted it's just missing from the signal. Here is an example where i wanted to upsample F=2 a cosine function:

The original signal is the plotted blue one and the upsampled is the red stemmed one. The last part of the signal is lost. The same goes for the downsampling by let's say F=5:

The blue one is again the original signal. In the first picture i needed to stretch the original so as to have a proper comparison and the same for the downsampled in the second picture. How can i fix this? Or is it fixable at all?


texts: Can I take piecewise FFT of a larger input sample?

texts:
I have a N point sequence x. Let X be the fft(x). (Assume x is arbitrary)

where,

x=[x1 x2 x3 x4 x5 x6 x7 x8]

fft(x)=[X1 X2 X3 X4 X5 X6 X7 X8]

How do I get the result of a 8 point DFT by doing cascading 4 point DFT and 2 point DFT?

Suppose I reshape the signal 'x' (of dimensions '1 x 8') to 'x1' (of dimensions '2 x 4')

x1=[x1 x2 x3 x4;
   x5 x6 x7 x8]

fft(x1 (first row)) i.e. fft(x1 x2 x3 x4) say =  [Y1 Y2 Y3 Y4]
fft(x1 (second row)) i.e. fft(x4 x5 x6 x7)  say = [Y5 Y6 Y7 Y8]

Then again taking fft column wise

fft(Y1 Y5)  say = [Z1 Z2]
fft(Y2 Y6)  say = [Z3 Z4]
fft(Y3 Y7)  say = [Z5 Z6]
fft(Y4 Y8)  say = [Z7 Z8]

However this Z will not be equal to X. How do I proceed to get X by doing these smaller point FFTs?

Any help would be great! Thanks


texts: Upsampling Methods for Computed-Tomography

texts:
I have two sets of data of given Field of view, one of them only covers a subset of the FOV of the other. I therefore want to upsample the one with the larger FOV to combine it with the other one.

So the question is: In the context of X-ray computed tomography, where the x-rays projections are backprojected to obtain a 3D reconstruction of the object, what is the best upsampling method for the purposes mentioned above?

Or to put it differently, what are the features in an x-ray image that ideally we don't want to damage too much... is it the contrast, the sharpness of the edges, etc. keep in mind, that as in many application of digital imaging, anti aliasing filters, Gaussian blurrings, etc. are performed on the images before them being processed for reconstruction purposes.

I know that each upsampling method has its merits, so considering the nature of the x-ray reconstruction methods, is there any way of deciding on one upsampling method based on the a-priori information of the object.


texts:
I have doubt, please clear.

A transform is said to be unitary transform when the product of the transform matrix and its conjugate transpose (Hermitian) matrix is equal to Identity matrix.

DFT for N= 4 gives

c =

 4     0     0     0
 0     4     0     0
 0     0     4     0
 0     0     0     4

Thus, it satisfies the Unitary property

But when N=2, we have

2      1+j
1-j     2

This doesn't satisfy unitary property right? Then how can we say the DFT is a Unitary Transform? Where I have misunderstood the concept? Help


texts: How do you perform range compression of SAR data on a 2D matrix?

texts:
I am building a Matlab program to process RADARSAT-1 data that came with this book.
The book describes the steps to creating SAR image as follows

I am currently implement the range compression step. For my matched filter computation I am performing a range FFT by transposing each row of the 2D matrix that represents my image first then taking the FFT of that.I understand what I need to do to perform my matched filter computation. My question is more related to how to do it in matlab on a 2D matrix. My 2D matrix  is 1536×2048 where 2048 is the number of columns which represents my fast time samples or range cells.Currently I am transposing each row so I can perform the fft on the rows of the 2D matrix in order to do the matched filter computation. Is that the right thing to do? If not what should I do?
Lastly, how can I confirm that my range compression code works?


texts: Matlab: Implement filter with no denominator

texts:
A filter like this:
$z^{4}\cdot\frac{2 + z^{-1} - 3z^{-3} + ...}{1 - 0.9z^{-1}}$ seems hard for me to implement. The fraction part is easy:

filtered_signal = filter([2 1 -3 ...], [1 -0.94], signal);

But how is the $z^{4}$ part being implemented?


texts: Speech synthesis requiring very little CPU performance?

texts: Back in the days of 1 MHz 8-bit CPU personal computers (Apple II, Atari 800, et.al.), there were software programs that could do comprehendible arbitrary text-to-speech synthesis on those PCs.  What published speech synthesis algorithms might be suitable for implementation on micro-controllers with similarly limited performance and memory?  (If any...)


texts: How to perform a Rubberband-Correction on spectroscopic data?

texts:
I basicaly understood how the Rubberband/Baseline-correction works.

The given spectrum is divided into (N) ranges.
The lowest points in every range are determined.
The initial baseline is built out of those points.
Now all the points on the spectrum are drawn down by the difference between the lowest point in the current range and the lowest point on the baseline.

There are some nuances, though, that I do not know how to handle. E.g., what if one of the points is exactly on the border between two ranges, etc.

Plus, I have to be able to prove that the algorithm that I am writing is a solid one and can be referenced by other works or scientific papers.

If anyone could give me some reference I would be very pleased.


texts: Closest two frequencies that can be resolved

texts:
Forgive me if I misuse the terminology.

I wish to transmit two slowly changing analogue signals by amplitude modulating two higher frequency sinusoidal carrier waves. I will superimpose the two AM signals, send them down a wire, and then demodulate them at the other end using a digital signal processor.

To demodulate them, I will simply use two DFT bins.

I want to select two carrier frequencies as close together as possible. What is the closest I can get two carrier frequencies without them interfering?


texts: Is Entropy Coding in Wavelet equivalent to Huffman Coding in DCT

texts:
How exactly we should describe the process of Entropy Coding in wavelets ? can't it be simply called Coding ?

Here is the block Diagram of Image compression using wavelet,

I am confused about the "Entropy Coding" part. Can i simply use Huffman coding to code bits ?


texts: How can I detect buttons in a GUI of an app with pattern detection using python?

texts:
I will explain taking as example the GUI of the Calculator app in windows.
Given this Calculator application open and focused, I need to find a way of detecting all the buttons.I can only use non-intrusive methods, so things like button id are out of the question. This leads me to image recognition. Or better say image detection, because I don't want a specific image, but a set of images that follow a certain pattern. I know how I can click/right-click/dbl-click/etc using a given button image and image recognition [1]. I don't know how I can do the other way round: scan the GUI and find those areas that meet the requirements of being a button (rectangle, text/icon/graphic labeled etc.). A bigger tackle would be to detect items that don't have the rectangle shape (e.g. icons on a windows desktop)

The nearest thing to what I need is detecting faces in a picture.[2] But I don't know how to apply this in my case. For human face detection, I saw that hundreds of pictures of the face must be used in the Haar cascade generation (don't know how I would do this with only 10-15 button snapshots). If another kind of object has to be detected like an apple, you need to generate Haar cascade for that object again using a lot of images.

Did any of you ever tried to detect buttons, items or whatever in a GUI using only pattern detection? I just need something to tell me "this is an icon/button" so I can grab that region in a snapshot.

[1] I use SikuliX with Python to perform actions on given patterns.

[2] I saw that this is done easily using OpenCV and Haar cascades (in XML format). Creating the Haar cacades requires a little bit of patience and skills.


texts:
Assuming that $f_{X|D}(x|d)=e^{\alpha d x - \psi(d)}f_0(x)$ for some constant $\alpha$ such that it is a valid pdf for every value of $d$, i want to establish that:
     $$\mathbb{E}[D|X=x]=\frac{1}{\alpha}\frac{\partial}{\partial x}\lambda(x)\quad \text{, where}\quad  \lambda(x)=\log\frac{f_X(x)}{f_0(x)}$$

*There's no missing information. It's from a collection of Information Theory problems.


texts: An intuitive explanation as to why a matched filter is time reversed?

texts:
It easy enough to study correlation, and matched filters. But the challenge I see unmet anywhere to date and which I struggle to meet myself is a simpler presentation, to with more difficult task I suspect.
Can you explain, in lay terms, to satisfy the intuitions of a listener, why the matched filter is a time reversed copy of the signal we hope to detect? It is not a mathematical proof, or explanation I'm after (they abound) but other some way of explaining it (to lay listeners) that generates a "a ha" experience, a feeling that it makes sense, and I get why we time reverse the signal for a matched filter ...
This is something I have yet to find, read or master myself.


texts: Is the calculation of variance correct from an estimator -- confusion regarding complex number where I am going wrong?

texts:
I have an expression for the varaince of measurement noise obtained from an estimator. The measurement noise is additive white gaussian having values in complex domain. I have squared the term $(.)^2$, should I be squaring or just implement without squaring? But if I dont square the term, then I am getting a complex valued variance which is incorrect as variance is always real valued.
Kindly suggest what to do.


texts: Sample text to collect all possible English biphones for Text-To-Speech

texts:
Since the list of phonemes in English is fixed, it should be possible to come up with a sample text(s) to collect all possible biphones for text-to-speech synthesis.

Does anyone have a sample text for this purpose?


texts: Why does signal averaging reduces noise levels by more than $\sqrt{n}$?

texts:
I have an electrophysiology signal, which is time-based. It basically measures neural activity in the form of potential differences over time.

The noise in this signal is assumed to be a random background noise caused by the device electronics and is assumed to have a normal (i.e., Gaussian) distribution (based on previous literature) and it is supposed to be stationary. Below I have added some plots that show a raw sweep (2.3 ms of recording per sweep, sample rate= 56k, 127 points were stored), the histogram of the values, and the FFT.

In answer to the answer below: I am fairly confident we can assume normally distributed data, since it's a bell shape. However, The FFT is not flat as answerer suggests it should be when normally distributed. However, we only have a few ms of data, which may complicate this type of analysis.

Assuming a normal distribution, theoretically, doubling the number of sweeps and averaging them should reduce the noise by a factor of the square root of $2$ (about $~1.41$). This kind of signal averaging to reduce random background noise is called ensemble averaging.

Now, I have a recording (without a signal, just background) and I have determined the noise level at $1$, $2$, $4$, $8$ and $16$ averages.
Edit: This was done by re-recording the background noise and determining the average with MATLAB, in the format B = mean(A1,A2), with  Ai being a vector with a time-based signal.

Noise level was defined as the standard deviation of the sweep. When I determine the factor of improvement between these averages, i.e., 2 vs 1, 4 vs 2, 8 vs 4, and 16 vs 8, I find the factors $1.4$, $2.5$, $2.5$ and $2.3$, respectively. These improvement factors were simply calculated by determining the ratio between SD$n$ / SD$2n$, with $n$ being the number of averaged sweeps.

The factors I find are averages across 31 electrodes measured in a total of N = 13 people. In other words, it's consistent across trials and 13 different devices. I didn't do stats, but I bet the latter three improvements are significantly $> 1.41$.

Why is the noise reduction higher than the theoretical 1.41 when doubling my sweeps from 2 to 4, from 4 to 8, and 8 to 16?

Raw sweep; 2.3 ms, 56 kHz sample rate, 127 data points

Histogram of amplitude distribution

FFT of the signal


texts: What is the x and y axis showing for this matplotlib.pyplot histogram of a .wav file?

texts:
The python code below generates a histogram from a .wav file, but what exactly does the histogram show? Is the x axis amplitude binned by sampling rate? Is the y axis a count of how many samples are in each amplitude bin? And how is amplitude calculated?
import numpy as np
import matplotlib.pyplot as plt

Fs, data = wavfile.read('audio file')

plt.hist(data, bins='auto')
plt.show()


texts: Relation between power spectral density and mean absolute value

texts:
The root mean square

$$\sigma_{x} = \sqrt{\frac{1}{T}\int_0^T x^2(t) \, \mathrm{d}t}$$

of a finite zero-mean random signal $x(t)$ in the range $0 < t < T$

is related to the signal's power spectral density $S_{xx}(\omega)$ via the relationship

$$\int_{-\infty}^{\infty} S_{xx}(\omega) \, \mathrm{d}\omega = \sigma_x^2$$

However, I was wondering if a similar relationship exists between the power spectral density and the mean of the absolute value of the signal, i.e.

$$\mathrm{mean}\left({|x|}\right) = \frac{1}{T}\int_0^T |x(t)| \, \mathrm{d}t$$


texts:
Context

I am reading this book Discrete Cosine Transform [...]. K.R. Rao and P.YIP. 1990.

Given That

From Equation (2.4.11.a)

$$
X^{C(2)}(m) = \left ( \frac{2}{N} \right )^{1/2} k_m \sum^{N-1}_{n=0} x(n) \cos \left [ \frac{(2n+1)m \pi}{2N} \right ]
$$

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

Here, we have

$$
k_p =
\begin{cases}
\frac{1}{\sqrt{2}} & \text{when } p = 0 \text{ or } N \\
1 & \text{otherwise.}
\end{cases}
$$

From Book's Equation (2.8.7)

$$
F^{C(2)}(k) = k_k \exp \left [ \frac{-j \pi k}{2N} \right ] \hat{F}_F(k)
$$

Here, we have

$$
k_p =
\begin{cases}
\frac{1}{\sqrt{2}} & p = 0 \\
1 & p = 1, ..., N-1 \\
0 & \text{otherwise} \\
\end{cases}
$$

Now Comes the following

  It is not difficult to show that the sequence $f(n)$ can be recovered from the transformed sequence $F^{C(2)}(k)$, as as a result we have

From (Equation 2.8.8)

$$
f(n) = \sqrt{\frac{2}{n}} \Re \left \{ \sum^{2N-1}_{k=0} k_k F^{C(2)} \exp \left [ \frac{j \pi k}{2N} \right ] \exp \left [ \frac{2j \pi k n}{2N} \right ] \right \}
$$

$n=0, ..., N-1$

Question

I do not understand how that $k_k$ is present on the (2.8.8). If I had to deduct this equation myself from (2.8.7) $k_k$ would pop on the other side of the equality inverted, as $k_k$, though it would yield a division by zero. But I could pretend (2.8.7) uses the same definition for as (2.4.11.a) resulting on a $F^{C(2)}(k)$ that is defined for $k > N-1$, then $k_k^{-1}$ wouldn't yield division by zero. How is (2.8.8) actually obtained?
Why having two different definitions for $k_k$?


texts:
I'm trying to understand what is the "phase coded modulation" and how does it work.

As far as I know, Although the linear frequency modulated (LFM)
waveform continues to be the workhorse of modern radars, there is growing interest in the use of
phase-coded (PC) waveforms.  Phase coding can be used to reduce radio frequency interference (RFI)
between adjacent radars and make the waveform less vulnerable to exploitation, and more robust to
interference.

and I'm not a familiar with like these radio frequency engineering. So.

Is the Phase modulation same "phase coded modulation"?
I'm confused that because there are several types modulation way in
https://kr.mathworks.com/help/comm/ug/digital-modulation.html#fp45657

  Kind of Modulation
  Amplitude modulation
  Phase modulation
  Frequency modulation
  Continuous phase modulation
  Trellis-coded modulation

But I can't find anywhere about Phase coded modulation.

Would you help me to understand what is the "phase coded modulation" and how does it work with very very simple example for referencing?


texts:
I have to draw a root locus with respect to K and find for which values of K the system is stable, then draw conclusions on the principal shape of the step response.

$$
G_{0}(s)=\frac{K(s+2)}{s(s+1)(s+3)}
$$

Isn't the step response just $\displaystyle \frac{1}{s}$?

I can get the characteristic polynomial of $G_c$ which is $s(s + 1)(s + 3) + K(s + 2)=0$, and I know it is in the form of $Q + KP$ wher $Q$ stands for $G_0$ numerator and $P$ for $G_0$ denominator, so the function is supposed to move between the zeroes to the poles.
I checked the solution and I still don't understand how to do the asymptotes, how do they come up with the intersection with the real axis, how do they draw the shape in form of a bow and 2 bars on the real axis. The intersection with the imaginary axis part makes much more sense.

What's the reasoning behind the asymptotes and the root locus?
The solution:


texts: Help understand the results of hough transform

texts:
I'm trying to apply the Hough transform of the following image using MATLAB:

The result is:

How could I predict that there are 6 bright spots in the resulted image? I expected 8 spots since I observe 8 lines in the original image.


texts: How to determine range of HSV values of the image?

texts:
Is there any method to find range of hsv values of image (object) in rectangle inside big image? I mean, there is fixed size rectangle on my image(frame from camera) and i want to find range of hsv values of object which is inside this rectangle. Object is changeable. If i put red ball inside rectangle it should find range of red ball's hsv values, if i put my hand inside rectangle it should find range of my skin's hsv values. Is there any way to do this?

Edit:
I get part of image with ROI, and convert it to hsv format. But problem is behind red ball, I get wide variety of colours. So which range will I choose from this image?

Example: my frame:

Now I want to find best range of hsv values of image in rectangle box for tracking it. (Now it should be green and ranges should be h:53-90 s:74-147 v: 160-255)

If i can find this values, by using inRange method, i can convert it to binary image for tracking. Like this:


texts: Transfer Function Block Diagram Confirmation

texts:
Can someone confirm if this is the correct block diagram for the following transfer function?

The original equation provided was: y[n+1] = y[n] + 0.01x[n]

Which I rearranged into H(z) = Y(z)/X(z) = 0.01 / (z - 1) and then also rearranged into y(n) = y(n-1) + 0.01x(n-1)

And left me with the following block diagram,


texts: Discussion of simulation results of a signal sequence: amplitude

texts:
I need your help in understanding if the simulation result is correct.
I have simulated a random BPSK sequence that was upsampled and filtered with a raised cosine filter.
The length of the filter is 150 taps, sampling frequency 16, cut-off frequency 0.5 and transition bandwidth is 0.2.
The noise was simulated with EbN0 = 10 with the following  noise power:
SignalEnergy = (trapz(abs(filtered_signal_tx).^2))*(1/Fs);
Eb = SignalEnergy/(2*Nb);
N0 = Eb./(10.^(EbN0/10));
NoisePower = 2*N0*Fs;

The signal has a phase offset(shift): exp(1j*(2*pi*t+phase_offset));
As a result, I have got a signal, which has a small amplitude. It is less than 0.08...
I don't know if it is wrong or correct?!
What do you think?


texts:
I want to segment the image below (click for bigger image) into columns based on the natural gaps that exist in the columns. The columns should be perfectly vertical and need not be the same size, but should split each column by itself.

I have tried Hough Line Transform, Scipy Label, convolution filters, specifying manual cuts, etc. all with very limited success and nothing I would consider a "solution".

What approach should I take to segment these types of images into columns? Preferably a solution that I can code myself or exists in OpenCV or Python.


texts:
If symmetry conditions are met, FIR filters have a linear phase. This is not true for IIR filters.

However, for what applications is it bad to apply filters that do not have this property and what would be the negative effect?


texts: How is the ends of a bandpass filtered signal reliable after applying the cosine taper?

texts:
A taper function decaying to zero at ends such as a cosine window is routinely applied to the signal before bandpass filtering to reduce end effects. Obviously the window alters the waveform.

So, as for the filtered waveform, how is it reliable near the beginning and end?

how does the the decay ratio of the window (e.g. 5%, 10% ...) influence? It's selected arbitrarily or by some criteria?

In the filtering functions of matlab and scipy, some extrapolation is applied to the ends of the signal to match their slope. In such situations, is the windowing still necessary?


texts: Performance evaluation in terms of channel MSE and BER

texts:
The first figure on the left (in this paper) provides the channel NMSE performance and the the second figure from left provides detector's BER performance.

I am curious, when the channel MSE overlapped in the first figure, why there is a difference between the BER performance.

I think, that there should be no difference, Moreover, even if there is a gap between MSE the use of quantizer of quantizer can bridge the gap in BER performance.

Thanks


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

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

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

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

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

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


texts: Determine the Unit Impulse Response of ARMA Filter for 8 First Indices

texts:
I have this problem, I've checked for almost 3 hours in books and internet for the unit impulse response for n = 0 : 7 but nothing.

I've benn working on how to get h[t] but I end up on a dead end at:

[

From there, I don't know how to solve this problem, Someone could help me pls?


texts: White Gaussian Noise: a reality or fantasy?

texts:
White noise has infinite variance, whereas Gaussian random variable must have some finite-variance to exist. How do two disjoint phenomena co-exist?

Help will be highly appreciated.


texts: Where do you recommend getting "standard" images for testing algorithms?

texts:
I've noticed that there are a LOT of websites which claim to have a standardized repository of images (for image processing algorithms). However, I don't know which ones are actually trusted by the research community, and which are not. I know some file formats hint of untrustworthiness (e.g. never use a JPEG since its compression is lossy, and it is probably not the original image). Many of these repositories appear to have the same image, but their actual files do not match pixel-for-pixel.
So, what is your recommended repository for getting standard images for testing image-processing algorithms (e.g. images like cameraman, boat, airfield, peppers, pirate, etc.)


texts:
Couldn't decide whether the question belongs on Mathematics or here, but I guess I'll ask it here.

I've got a discrete - valued function defined on a finite segment of a line that looks like this:

Then the task is, given a set of other functions, find a function in this set that is the closest - in some way that I'm about to describe - to the original one.

The problem is that the candidate functions may be, actually, circularily shifted with respect to one another. If the functions were not shifted, any reasonable norm of difference (l1, l2, any lp, I guess) would've been okay. I also want to notice that I do not necessarily want a function that satisfies the mathematical definition of a norm (or even a distance) - I just want the problem solved.

The obvious thing to do is to look at the maximum of circular cross - correlation. However, this is actually a bad idea because even if I normalize both compared function by their max values, the function f(x) = 1 will give the maximal value with any other function.

I want to know whether there is some standard approach for such a situation or should I just compute an l1 - norm for every possible circular shift and then take the minimum value.


texts: Algorithm for testing input sound for multiple expected notes?

texts:
** originally I asked this question on stack exchange and was pointed to dsp.stackexchange **

I need to determine when the correct musical notes (from a piano) such as in a chord are being played. While detecting single notes is a no brainer, and full music transcription seems to still be considered a research problem, I'm hoping that since my program knows which combination of notes to expect and only needs to classify the played notes (maximum of six) as correct/incorrect that the task will fall somewhere in between.

Can anyone point me towards an efficient algorithm (I need to implement it in JavaScript) that should be able to achieve this with minimal false negatives (incorrectly telling the user they hit the wrong note)?

My current idea is to test that the signal has better correlation with the target frequency than the nearest semitone above and below for each expected note.


texts:
I'm trying to create a basic digital signal processing pipeline, to manipulate a WAV file.
Nothing fancy, I'm just looking to apply some filtering.

My first attempt was to divide the input signal to N rectangular windows (which do not overlap), convert each window to the frequency domain, apply the filtering of my choice to the spectrum, convert it back to the time domain and insert the resulting window in place of the original one.

(By rectangular window, I mean rectangular in the time domain.)

This worked, however, the only downside was that there was a discontinuity at the edges of each window, severe enough to distort the overall signal and introduce harmonics of significant amplitude.

My attempt to resolve this was to replace the rectangular windows with Hann windows, and apply the same process.
With non-overlapping windows, the obvious result was that the output signal became a set of Hann "pulses", which is as bad as the discontinuities.

I seem to be missing a very basic part of digital signal processing.

I understand the windows should overlap, but how exactly? Should the overlapping portions be averaged? How do we ensure that the amplitude of the resulting signal is unchanged despite Hann window overlapping/averaging?


texts: Issue with the time vector returned by $\tt signal.spectrogram$ function

texts:
I'm puzzled by some very simple concept with building up a spectrogram.

Here is a toy example of the issue:

import numpy as np
import scipy
from scipy import signal
import matplotlib.pyplot as plt

f0 = 5
t = 2
fs = 250
N = t*fs
t = np.linspace(0, t, N, endpoint=False)
x = np.sin(2*scipy.pi*f0*t)
f, t, Sxx = signal.spectrogram(x, fs)
plt.pcolormesh(t, f, Sxx)
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [sec]')
plt.show()

What I get for t puzzles me, I get [0.512,1.408]. I would expect to get an array starting from 0.0 to 2.0 seconds long - so if I want to see the power spectrum density in 1.8 seconds I could easily do so.

I probably miss something pretty basic, so it would be great if someone could shed some light on this.

Mike.


texts: I generate a 19kHz audio signal with sampling frequency of 44100 using Matlab. Due to the sampling, some other frequencies are generated in addition to the 19 kHz signal. When I play this audio, I hear additional voices. How I can remove these additional sides? I think it needs aliasing filtering, but I cannot filter signal after DAC, because I generate the audio signal and play it with phone or computer, I cannot put any filter after A/D.


texts: What should I know about Morphological processing

texts:
i'm studying my first year as technical student , and we are studying image processing, so last week we got delivered a topic to work on : Morphological processing /opening
, i'm trying to dig deeper in internet but i lost in the amount of informations as a beginner, I still can't figure out what is that (morphology processing) all about.
Can you give me a hand with this?


texts: How to check if a signal is power signal or energy signal?

texts:
What is the procedure to check the signal type?

example:

$ x(t) = A  \sin (\omega t) $

$ y(t) = A e^ {-\lambda |t|} $


texts:
I have a discrete signal that goes as follows:
$$x[n]=[-1,4,-1,0]$$
I have already done the DFT for the signal, with the following result:
$$X[0] = 2, X[1]=-4i,X[2]=-6,X[3]=4i$$
But for some reason, I seem to have trouble getting the frequency resolution for these samples(?).
My sampling interval (T) is 0.03125 seconds, and with some logical deducting I came up with these following frequencies:
$$f_0=0,f_1=8,f_2=16,f_3=8$$
However, I do not understand how I came to this conclusion, but it seems correct. I tried to use the following formula:
$$\Delta f = \frac{1}{T_{0}} = \frac{1}{NT}$$
Could someone explain how I ended up with these frequencies?


texts: How did they color old black and white movies?

texts:
This question is on dsp.SE as I'm mostly interested in the signal processing part.

There is an Indian movie Mughal-e-Azam which was released in 1960 in black & white which has been reproduced in color in 2004.

How did they color each pixel perfectly?
What technique have they
used to identify the color placement on each pixel?

Look at one of the screenshots from the movie:

I've got an Einstein black-and-white photo which I want to colorize. How is it possible to do so  without knowing what he was wearing back then and what was the actual color of his clothes, background etc.


texts: How to define a Mel Spectrogram bank in terms of frequencies?

texts:
I am trying to determine the formulation of a single Mel Spectrogram bank in terms of frequencies. For ex. in the following image I have 256x80 banks in total, and I would like to determine what would be the mathematical description of a given column i.e. 1x80 banks.

In the following plot I have set everything to 0 and am varying frequency 1Hz from -1 to 5. It is emperically evident that frequencies and banks have many-to-many mapping and one bank is a function of 0 or more frequncies, and any given frequency appears in 0 or more banks.

I believe the formulation might be something like:
$$Mel_{j} = \sum_{i}^{N}{f_{ij}(F_{i})}$$
Where $Mel_{j}$ is one of the mel banks, $F_{i}$ is the value of frequency $i$ and $f_{ij}$ is a function that describes mapping of $F_{i}$ to $Mel_{j}$
I was wondering if this formulation is well known and if not someone might have worked on something similar and could help me figure it out.


texts:
I am trying to use FFTW function for IFFT:

fftw_plan_dft_c2r_1d(size, out_cpx, out, FFTW_ESTIMATE).

My goal is to input desired frequency response data and obtain related coefficients using Inverse Fourier Transform. Since out_cpx is complex, I calculated magnitude (0,1) and phase, taking into account it is linear phase filter. Since I know the phase I calculated real and imaginary parts as Acos(w) and Asin(w) and with these pairs I populated out_cpx[i][0] and out_cpx[i][1]. When I perform Inverse Fourier Transform using fftw_plan_dft_c2r_1d() I assume I am right away getting FIR coefficients h[n] as the result. With these coefficients i go to MATLAB and get response using freqz(C,1) where C holds my coefficients. But I am getting a lot of ripples in my graph for response and for phase. Am I doing something wrong here? It is not clear to me how to specify phase of the desired frequency response.


texts:
I want to implement the 'filter' function in matlab but I just can't seem to replicate the results I get when using the matlab function.

My understanding of the matlab function is that it takes 3 arguments: 'b' = array of numerator coefficients , 'a' = array of denominator coefficients and 'array' = the array to be filtered.So if you wanted to filter an array with a filter H(z) where:

$$
  H(z) = \frac{b_1 + b_2z^{-1}}{a_1 + a_2z^{-1}}
$$

and store it in an array called "result" then you would write

result = filter( [b1 b2] , [a1 a2] , array );

Now I tried to implement the 'filter' function ( with the number of coefficients in the numerator and denominator limited to 2) by getting the difference equation from H(z) then using it to get the resultant array but this route does not seem to work, ie

$$
  H(z) = \frac{Y(z)}{X(z)} = \frac{b_1 + b_2z^{-1}}{a_1 + a_2z^{-1}}
$$

$$
  X(z) (b_1 + b_2z^{-1}) = Y(z)(a_1 + a_2z^{-1})
$$

$$
b_1x[n] + b_2x[n-1] = a_1y[n] + a_2y[n-1]
$$

$$
\implies y[n] = \frac{b_1}{a_1}x[n] + \frac{b_2}{a_1}x[n-1] - \frac{a_2}{a_1}y[n-1]
$$

Meaning the matlab code should look something like

function [ result ] = myFilter( b , a , audio )
%MyFilter
%a,b = Two Element Arrays
%array = the array we want to filter

N = size(audio,2);

for i = 1:N
   if( i == 1 )
       result(i) = ( b(2)/a(1)) * audio(i);

   else
       result(i) = ( b(1) / a(1) ) * audio(i) + ...
                   ( b(2) / a(1) ) * audio(i-1) - ...
                   ( a(2) / a(1) ) * result(i-1);
   end
end

end

The function above does not produce the same results as the matlab filter function.Does anyone know what I am doing wrong? Have I understood how the filter function in matlab works?


texts: Custom ringtone sounds good on computer, can't be heard on phone

texts:
I did a ringtone using Audacity and then used this website to convert to mp3. I sent this mp3 file to a friend who said he couldn't hear it well. I have an iPhone so I converted it to m4r. Both the mp3 and m4r files play nicely in my computer. However, when I imported it to my iPhone SE, only certain parts could be heard. Here it is if you want to try it

https://soundcloud.com/jpproducoes/test-ringtone/s-TmLZ7

Open it on your computer, then open the same link on your phone (you can use soundcloud, the result will be the same). On the phone it will sound much more distorted, and some parts won't be audible. I believe this might be due to sampling rate problems, but I am not sure. Any ideas?


texts: Difference between slow, medium, fast AGC

texts:
In digital receivers, at the RF level (not audio), I often read about "slow AGC" and "fast AGC", and sometimes even "medium AGC" (but I think the latter is for audio). From my understanding, the slow AGC reacts to changes in noise level and variations in signal level over a relatively long period of time (fading for example). And from what I understand, the fast AGC is used more on symbol level, trying to keep them at a relatively stable power level, but not too fast so that amplitude-dependent symbols would not be affected (like M-QAM or M-PSK with M > 4).

Is there more to it? Also, would the "medium AGC" have any use in the RF chain?


texts: Residual error when setting measurement matrix in compresssive sensing

texts:
I have an issue when implementing compressive sensing to recover sparse vector. Assume I have sparse vector $x$ of length, for example, $(256,1)$. $x = [x_1,x_2,.....x_{256}]$. This vector is transformed into time domain using Fourier matrix, $X = iFFT(x)$, and then convoluted with a channel $h$ resulting $y = X*h$, where * denotes the convolution operation. Let's make $H$ denote the toeplitz matrix of vector $h$ correspondent to $y = XH$

$case-1: $ Let's describe the above process as below:

$1). $ $X = F'x$ , where $F$ is Fourier matrix transformation obtained in matlab by F = dftmtx(256). and then $F'$ is the inverse Fourier matrix. This $X$ is correspondent to X = ifft(x).

$2).$ The signal $y = HX$ can be expressed in function of $x$ as $y = HF'x$. The use of compressive sensing here is straightforward by putting $y = Ax$ where $A= HF'$ is the measurement matrix. Therefore, sparse vector $x$ can be recovered using any algorithm such as OMP, MP or any other.

$case-2: $ Assume that we transformed the vector $x$ into time domain using $iFFT$ operation, but his time in vector of $32$. It means that we take every $32$ elements of $x$ and then perform $iFFT$ for it. This process can be performed in matlab as, X1 = reshape(x,32,[]); X_2 = reshape(ifft(x),[],1); This $X_2$ is equivalent into the $X$ mentioned in case_1  but with Fourier transformation of vector $32$

$1)$ Similarly to above manipulation, $X_2 = F'_{32} x$, where $F'_{32}$ can be expressed as below :

where $X_2$ here is correspondent into $X$ in the case_1, and $F'_{32}$ is inverse Fourier matrix of size $32$ obtained in matlab as dftmtx(32)' .

$2)$ The signal $y$ can be expressed similarly to case_1 as $y = HX_2 = HF'_{32} x$. BUT, performing compressive sensing into this case, $y = A_2x$  where $A_2 =  HF'_{32}$ doesn't recover the vector $x$.

My question, theoretically, any compressive sensing algorithm such as OMP, MP should be able to recover $x$ in both case since measurement matrices are known and then should be straightforward, but in case $2$ where measurement matrix $A_2$, the CS doesn't work? A small residual error happens when performing that in matlab, so is there mathematical expression for that error? Is it possible to overcome that issue ?

EDIT:
Here is the method how to create the $A_1$ and $A_2$ in matlab

F = dftmtx(256)/256;     %%Fourier matrix
X = F*x;                 %%x is sparse vector of length(256,1) and X time domain of x
h = randn(32,1);         %%channel
H = toeplitz([h(1) zeros(1,255) ], [h.' zeros(1,255) ]).'  %%toeplitx matrix
y = H*X;                 %%recieved signale

A1 = H*F;               %measurement matrix
x_hat = omp(A1,y,q)    %estimated x, and q is the number of sparse

%%% creating A2
F2 = kron(eye(8), dftmtx(32)/32);   %Fourier matrix
X2 = F2*x;                % ifft of reshaped x
y2 = H*X2;

A2 = H*F2;
x_hat = omp(A2,y2,q)    %estimated x, and q is the number of sparse

where omp is copressive sensing OMP algorithm. In case of building the measurement matrix using $A_1$, x_hat can be recovered well but when using the $A_2$, a small error occurs, it means that estimated signal can't be demodulated and have BER = 0 compared with original x

Here you can get the complete code Here
Thank you


texts: IIR: get magnitude and phase from coefficients

texts:
i'm trying to plot 1)the magnitude and 2)the phase response of an IIR eq with only its coefficients.

i've found somewhere that one can find the magnitude response starting with the transfer function :

  ex: y(n) = a1*y(n-1) + b0*x(n) + b1*x(n-1) + b2*x(n-2)

  ... has the transfer function

       Y(Z)    b0 + b1 Z^-1 + b2 Z^-2
H(Z) = ---- = ----------------------
       X(Z)     1 - a1 Z^-1

  To find the frequency response at theta=pi (folding frequency)

  Recall Z = exp( j*theta )

  Z = exp(j*pi) = -1

         b0 + b1*(-1) + b2*1
H(-1) = ----------------------
            1 - a1*(-1)

let's try with a frequency of 0.25*Pi

exp(0.25*j*pi) = 1 + j*sin(0.25*pi) ~= 1 + 0.72*j

but after some calculations i end up with a complex division:

     (b0 + b1 + 0,52*b2) + j*0.72*b1
  =  -------------------------------------
      (1 - a1) + a1*0.72*j

so where do i go from there?

and how to compute the phase ?

thanks a  lot!

Jeff


texts: Determine if $ y[n] = ny[n-1] + x[n]$ is linear time invariant and BIBO stable

texts:
Check if the following system is linear time invariant and BIBO stable..
$$
y[n] = ny[n-1] + x[n]
$$

for $n\ge 0$. We are also given that the system is at rest (i.e. $y[−1] = 0$).

I know that to show linearity, we show that
$
H[ax_1 + bx_2] = aH[x_1] + bH[x_2]
$
except that the question has a $y[n-1]$ on the RHS which is throwing me off.

Edit: My solution:
$
LHS = $ $$H[ax_1[n] + bx_2[n]]
$$$$
= n(ay_1[n−1] + by_2[n−1]) + (ax_1[n]+bx_2[n])
$$$
 RHS =$$$  aH[x_1] + bH[x_2]
$$$$
=  any_1[n−1]+ax_1[n]+bny_2[n−1]+bx_2[n]
$$
So $
 RHS = LHS
$

Edit 2: New issue when I try and show time invariance.

If $ x_2[n] = x1[n-k] $ then prove $y_2[n]  = y_1[n-k]$.

BUT $RHS= $
$$
(n-k)y_1[n-k-1] + x_1[n-k]
$$ And $LHS = $
$$ny_2[n-1] + x_2[n] = ny_2[n-1] + x_1[n-k]$$


texts: Is there a hypothetical situation where once can use Wiener filter prior Matched filter to get better result? Mostly I am looking for an example or a situation where both the filters are required as ther are doing different thing to the received signal. Under what circumstanctes this would be helpful?


texts: Understanding and converting Matlab's coarse frequency compensation function

texts:
I am trying to understand and convert to C Matlab's coarse frequency compensator the documentation for which can be found here: Documentation

Matlab's design itself is based on this IEEE paper: Carrier Frequency Recovery in All-Digital Modems for Burst-Mode Transmissions

I want to translate the correlation based algorithm for QPSK modulation; for example the invocation of such an object in Matlab would look like this:

cfc = comm.CoarseFrequencyCompensator('Modulation', 'QPSK', ...
                                      'Algorithm', 'Correlation-based', ...
                                      'MaximumFrequencyOffset', 12.5e3, ...
                                      'SampleRate',400000);

The documentation states that the estimate is given by:

$estimate=\frac{f_{samp}}{\pi(L+1)}arg\{\sum_{k=1}^{L}R(k)\}$

where

$L = round(\frac{f_{samp}}{f_{max}}) - 1$

$R(k) = \frac{1}{N-k}\sum_{i=k+1}^{N}r_{i}r^{*}_{i-k}$, autocorrelation function,

$r_k = e^{j(2\pi\Delta f k T_{s} + \theta)}$,  model of received signal

$r^{*}_{k} = e^{-j(2\pi\Delta f k T_{s} + \theta)}$, is the complex conjugate of $r_{k}$

$N = \frac{L+1}{M} $, number of samples

$M = 4$, because for QPSK there are four possible symbols

My C code:

double complex autocorrelation_func(double complex* input,
                                    double complex* conjinput,
                                    int k, int input_size)
{
    double complex auto_corr_sum = 0+0*I;

    for(int i = 0 + k; i < input_size; i++ )
    {
        auto_corr_sum += input[i]*conjinput[i];
    }
    auto_corr_sum *= (1/(input_size-k));

    return auto_corr_sum;
}

double coarse_freq_comp(double complex* input,
                        double complex* output,
                        int input_size, int mod_order,
                        double max_freq_offset,
                        double sample_rate )
{
    int auto_corr_seq = round(sample_rate/max_freq_offset)-1;
    int num_samples = (auto_corr_seq+1)/mod_order;
    double mult_const = sample_rate/((auto_corr_seq+1)*M_PI);
    double complex input_conjugated[input_size];
    double complex auto_corr_sum = 0 + 0*I;
    double offset_est = 0;

    for(int ndx = 0; ndx < input_size; ndx++)
    {
        input_conjugated[ndx] = conj(input[ndx]);
    }

    for(int i = 0; i < num_samples; i++)
    {
        auto_corr_sum += autocorrelation_func(input, input_conjugated, i, input_size);
    }

    offset_est = mult_const*catan(auto_corr_sum);

    return offset_est;
}

Questions:

How do I apply the estimate to the input signal to get a compensated output signal?

What is wrong with my C code? I am not getting the same estimate value that I get for a given complex array in Matlab. I must be doing something wrong but this isn't obvious to me.


texts: Autocorrelation - Understanding reduced correlation at periodic time shifts using np.correlate (versus statistical autocorrelation)?

texts:
I'm going through the Think DSP by Allen B. Downey, and I'm struggling to understand a specific aspect of np.correlate and how it differs from statistic autocorrelation. The question is at the bottom, but to give it here: why does the magnitude of autocorrelation drop off at periodic intervals of time shift using np.correlate but not the statistical autocorrelation function (pearson product)?

The presented signal is

Using statistical autocorrelation, which I believe normalizes the signal and unbiases it, it looks like:

Using np.correlate(signal, signal, "same") I get

The book mentions that "you’ll notice that the correlations computed by np.correlate get smaller as the lags increase. That’s because np.correlate uses the unstandardized definition of correlation; as the lag gets bigger, the number of points in the overlap between the two signals gets smaller, so the magnitude of the correlations decreases." - This is the part I do not understand. It's effectively a periodic signal, so why does the number of overlapping points get smaller with greater time shift? At some point the shift matches the period, and the correlation should be near 30 again? Does it have something to do with the lack of normalization in the latter?


texts: Signal equation using signal waveform and Fourier series?

texts:
Part A: Obtain signal waveform from mathematical equation of the signal

Let a sinusoidal periodical signal is represented by an equation

$$y=f(t)=10+10\cos\left(\frac{2\pi f_1t}{T} +\frac{\pi}{6}\right)+10\sin\left(\frac{2\pi f_2t}{T}+ \frac{\pi}{3}\right)\tag{1}$$.

Here, let us take $f_1=10 ,f_2=5 ,T=100$

Now, with this we can get the waveform of the signal $y$ on Cartesian axis with magnitude as $Y$-axis and time as $X$-axis.

(Sorry but I don't have any software tool now so that I could show the waveform diagram)

Part B: Obtain mathematical equation of the signal from its waveform and Fourier series formula.

i.e. part B is reverse of part A.

Here, my question is that how can we obtain (probably the same) mathematical equation $(1)$ for signal $y$ from waveform of signal $y$ and Fourier series formula
?

Fourier series formula for any sinusoidal periodical signal $y$ is given as

$$ y(t) = \sum\limits_{m=0}^{+\infty} a_m \cos(w_m t) + \sum\limits_{m=0}^{+\infty}b_m \sin(w_m t). \tag{2}$$


texts: When Is a Kalman Filter Different from a Moving Average?

texts:
this thread asks when a discrete time Kalman filter is better/different from a simple moving average of the observations:

Why use a Kalman filter instead of keeping a running average?

there's no definitive answer. can someone give a definitive example where the kalman filter, ideally in simple 1D case, does something different (and better) than keeping a moving average, and state the conditions when the kalman filter would reduce to a simple moving average?

one thought is that the kalman filter would not weigh all the data points equally because its variance is initially smaller and gets better with time. but it sounds like that would only matter close to initial observations and that once the variance converged, the kalman filter would weigh each observation equally just like a moving average, so don't see when the two are different and why/when the filter would do better.


texts: Combine FFT's of shorter length than sample data to get spectrum of all data

texts:
I have a set of sample data, a logarithmic swept sine as it happens, that is 1 second in duration sampled at 100kHz giving 100k samples. I want to calculate the frequency spectrum of the complete data however the FFT library function available only supports 4096 samples maximum.

Is it valid to sequentially FFT 4096 samples and sum the resulting data. Thereby ending up with 4096/2 frequency points that show the frequency content of the entire 100k samples?

I guess it is necessary to window each block of 4096 samples to be FFT'd but this would mean that the portions of the frequency content would be reduce in amplitude artificially.

How does one handle the phase information?


texts: Notch filtering minimizing artifacts in signal

texts:
I have to notch filter the following signal

As you can see it is a signal with really high amplitude.And this causes me a lot of problems. I've designed a bandstop filter (46-55 Hz).The MATLAB code is

Fs=5000; sampling frequency

Wp_NOTCH = [36/(Fs/2) 71/(Fs/2)];
Ws_NOTCH= [46/(Fs/2) 55/(Fs/2)];
Rp_db=-20*log10(.95);
Rs_db=-20*log10(.05);
[order_NOTCH,wn_NOTCH] = buttord(Wp_NOTCH,Ws_NOTCH,Rp_db,Rs_db);
[B_NOTCH,A_NOTCH]=butter(order_NOTCH,wn_NOTCH,'stop');

And the filter is
When i filter the signal (using filtfilt) I've the following signal As you can see it has a lot of artifacts,and "waves". Is there a way to minimize these artifacts?How can I delete the 50Hz component minimizing these artifacts?Any suggestion? Thanks in advance


texts: Algorithm to estimate coherence when measuring a transfer function with a swept sine

texts: I use a swept sine to measure a transfer function, multiplying the reponse by the in phase and quadrature excitation to get a complex response at each frequency. Can you suggest an enhancement to this technique to estimate coherence as well to give an estimate of measurement quality, especially at higher frequencies where the response is small.


texts: I just read the topic about Ergodicity but I have ambiguity about its meaning (by intuition). What does mean: (for mean) Statistical average = Time average. Could you please explain it in detail. Thanks.


texts: FIR difference equation vs convolution summation

texts:
I'm sorry if my question is too simple, but I am self-studying FIR because we have to report it in our class. I am using the book of Proakis and Manolakis as reference.

Under 7.2 Structures for FIR systems:

  In general, an FIR system is described by the difference equation
    (7.2.1)

  or evidently, by the system function

   (7.2.2)

Under 7.2.1 Direct-Form Structure:

  The direct-form realization follows immediately from the nonrecursive difference
  equation given by (7.2.1) or, equivalently, by the convolution summation

   (7.2.4)

I can see that they almost look the same apart from the bk in 7.2.1 and h(k) in 7.2.4. Question is, how are they different or related to each other? Sorry if you find my question too basic. Thank you in advance for the help!


texts: Sampling signal after root raised cosine

texts:
I am generating series of 1 and 0. then pulse shaped them with root raised cosine and after than deciding to do match filtering and recovering the bits however it turns out that I am not getting the correct bits.
What I do is as follow:

generate 0 and 1.
pulse shape them with RRC.
match filter using the same transmited filter
Remove samples up to peak using [val,pos]=max(mysignal)
sample using mysingal(pos:overample:end)

But what I get is not correct. I simulated this using 10 bits but I get only 2. I know there is something wrong but I do not know why. My method works when I only have one bit.

I am wondering if I have more than one bit should I do this procedure for every  bit.How can I approach this problem.

This is also my code assume RX2 is the pulse shaped version of my signal and os is the oversampling factor and 0.25 is the roll of factor. I am wondering how to detect all the bits like all 10 or 100?

RX1=filter(srrc(os,0.25),1,RX2);
[val,pos]=max(RX1);
RX=RX1(pos:os:end);
Bits=RX;
One=Bits>=Th;
BER=sum(One)/length(Bits);

The first one is the signal after passing through matched filter using conv command and no built in command like filter. the original bits are 0 0 1 1 1
Can anyone tell me how can I sample signal to get these bits back.
what command I should use. simply downsampling will not work.

The second one is for the same bits but using filter command for 0 0 1 1 1 can anyone please tell me how to obtain correct bits back from any of these figures.


texts: Simple VAD (Voice activity Detector) implementation

texts:
I am to classify few recorded phone calls as either containing speech or not by sampling lets say first 10 seconds.

The files contain either silence, dialtone or ringtones, or real human voice.

I have tried implementing butter-worth filter 100hz-400hz, then calculating Short time energy, Zero crossing rate, then calculating variance of the resultant array. But the results aren't good, calls with ringtones are confused with human voice.

I am implementing this in matlab or python (scipy & numpy)


texts:
I would like to take a time series of several spikes or narrow Gaussians and bandpass using an ift/ifft together with a frequency domain mask to create a series containing band-limited spikes. This is pretty simple but I have an additional desire - that the resulting band-limited spikes are causal - that is each band limited spike has a sharp onset (at the position of the original spike) followed by a tail dependent on the bandpass parameters.

The ultimate goal is a series that looks superficially like a superposition of a multitude of earthquakes.

Is there a simple way to do this?

Thanks very much for any guidance.


texts: partial fractions expansion inverse Z-transform, help

texts:
I have the correct solution from teacher's solution guide, but I was slightly confused by some algebra about the partial fractions expansion evidently

difference equation is as follows

$ y[n] = \frac{1}{4}  \cdot \left( x[n] +x[n-1]+x[n-2]+x[n-3] \right) $

where $x[n]$ will be unit-step sequence $u[n]$

PROBLEM STATEMENT:

Find what will the output sequence y[n] when input x[n] = u[n] = unit step

We can use tactic to get X(z) and multiply X(z)*H(z) so we get Y(z)
Then inverse Z-transform to get y[k]

we know from earlier calculation that the transfer function for this system willbe as follows

$
H(z) = \frac{z^3 + z^2 + z +1}{4z^3} = \frac{1}{4} * \frac{z^3 + z^2 + z +1}{z^3}
$

because unit step input sequence

$
X(z) = \frac{z}{z-1}
$

We can find out  $Y(z) = X(z) * H(z)$

I had some difficulty with the partial fraction expansion because of the $\frac{1}{4}$ multiplier, how does that constant multiplier term become accounted into the partial fraction expansion?

We were tipped by our teacher to find out $Y(z)/z $  , so that we can find good matches in our z-transforms-table. So, I try to make  Y(z)/z into the partial fraction expansion...

$
Y(z)/z = 1/4 * \frac{z^3 + z^2 + z +1}{z^3*(z-1)}
$

So, do we just "ignore" the constant multiplier $\frac{1}{4}$ on the r.h.s of the partial fraction as follows???

$
Y(z)/z = 1/4 * \frac{z^3 + z^2 + z +1}{z^3*(z-1)} = \frac{A}{z} +  \frac{B}{z^2} +  \frac{C}{z^3} +  \frac{D}{z-1}
$

So, after the "equating coefficients technique" I obtain followins system of equations

4A + 4D = 1

-4A + 4B = 1

-4B +4C = 1

1 = -4C

from these I obtianed coefficients as follows

A= -3/2

B= -1/2

C= -1/4

D= 1

With these results we should multiply both sides by z, so we get Y(z) in partial fractions format

$
Y(z) = A + B/z + \frac{C}{z^2} + \frac{D*z}{z-1}
$

From here, we should be able to find inverse-Z-transforms from the table as follows

$
Z^{-1}(Y(z)) =  -3/4 * Z^{-1}(1) - 1/2 * Z^{-1}(1*z^{-1}) - 1/4 * Z^{-1}( 1 * z^{ -2 } ) + Z^{-1}( \frac{ z }{ z-1 } ))
$

$
y[k] = -3/4 * \delta[k] - 1/2 * \delta[k-1] - 1/4 * \delta[k-2] +u[k]
$

where $\delta[k]$, is the unit impulse, aka the deltafunction

I put that into my Texas instruments calculator and the function seemed to hold true for some initial values

expected values are that for negatives indexes output will be zero

then for k=0: output = 1/4

k=1: output = 2/4

k=2: output = 3/4

k>=3: output = 4/4


texts: Wavelet is a small wave, what is problem with using sin(x) of one cycle?

texts:
As far as I know a $\sin(x)$ which is one cycle long can be made into a wavelet and then we can use its dilated and translated versions to represent another function, effectively a wavelet transform.

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


texts: How to implement an RF Mixer for chirp signal in python?

texts:
According to this link, all I need to do to mix two signals is to multiply them.
I have two chirp signals with same properties, one leading the other. If I multiply them to get the IF output, I don't get a constant frequency which is what I expect, I get double the frequency:

from scipy.signal import chirp, spectrogram
import matplotlib.pyplot as plt
import numpy as np
import math

def plot_spectrogram(title, w, fs):
  ff, tt, Sxx = spectrogram(w, fs=fs)
  plt.pcolormesh(tt, ff, Sxx, cmap='gray_r')
  plt.title(title)
  plt.xlabel('Second(s)')
  plt.ylabel('Hz')
  plt.grid()

fig, ax = plt.subplots()

T = 0.4
fs = 7200
t = np.arange(0, int(T*fs)) / fs
f0, f1 = 5, 120
sc = chirp(t, f0=f0, f1=f1, t1=T)
rc = chirp(t, f0=f0, f1=f1, t1=T, phi=90.)

plt.figure(1)
ax.plot(t, sc, 'blue', lw=1.25, label='sent')
ax.plot(t, rc, 'green', lw=1.25, label='recv')

plt.figure(2)
plot_spectrogram(f'Chirp mul, f(0)={f0}, f({T})={f1}', sc * rc, fs)

plt.show()


texts:
Is the following equation established?

  $$\int_{-\infty}^{\infty}s(t)r^*(t)dt=\int_W S(f)R^*(f)df $$

  $$s(t)\xrightarrow{\text{Fourier Transform}}S(f) $$

  $$r(t)\xrightarrow{\text{Fourier Transform}}R(f) $$


texts:
I'm trying to evaluate received signals from a transceiver. I want to make a threshold of received signals so that we can say some signals are meaningless.

I have learnt something about Signal-to-Noise-Ratio (SNR) and Received-Signal-Strength (RSS). But I found that symbols can still be decoded even in low SNR ($\le 0$ dB).

Here are my questions:

If I choose SNR as my metrics, will my threshold be relevant to modulation or coding schemes?
Will RSS be suitable for theoretic analysis as metrics?


texts: calculate and compare Ballistocardiogram (BCG) with Photoplethysmogram (PPG)

texts:
I have a task from my boss on calculating Ballistocardiogram (BCG) with Photoplethysmogram (PPG) using MATLAB. I would like to know how it can be done. Any kind soul can guide me?

As for what I know about how the test will be like is that I will need to take the BCG data and compare it against the PPG (I know Matlab has comparing function). So, in short, the PPG will be like a benchmark.

Ideally for BCG and PPG should have the same amount of data, but let's say the BCG for example 12 data and PPG only have 10. How do we compare it? because at every point it should match and I can't be removing the 2 extra data as it would be manipulating and that defeats the purpose already.

BCG: Ballistocardiogram,
PPG: Photoplethysmogram

I hope you guys can show me some light

-I don't have enough reputation to upload pic-

https://skydrive.live.com/embed?cid=0525DA685954952E&resid=525DA685954952E%21405&authkey=AOrI6qD7kFWV4_g
Diagram 1: find the differences of the corresponding data. As you can see in this diagram 1B against 2B

https://skydrive.live.com/embed?cid=0525DA685954952E&resid=525DA685954952E%21406&authkey=ACIhYwrTqZaNnp8
Diagram 2: let's say BCG has 12 data and PPG only 10(hopefully this don't occur). In this diagram, should 1D correspond to 2C? or? and if we remove 1C this will be manipulating the data.

if there are other suggestions please tell me also. I don't need to use Matlab to compare. As long as I can compare, anyway is ok


texts: Carrier Recovery, Frequency and Phase offset in (CP)FSK

texts:
I wanted to ask what exactly is meant for phase recovery of signal (at least intuitively). I have a working example of an acoustic fsk modulation and demodulation software but i don't know how to apply these techniques to improve the detection of the signal.  The center frequency is 30Khz and demodulation is done in passband.

In the following i will briefly explain how it works:

1- Send a chirp, used from the receiver to find the beginning of the signal

2- Wait for a fixed amount of time, known also to the receiver

3- Send the modulated  signal

4- Send a final chirp

The receiver after i get the  correlation v with the chirp, i buffer data, isolate the signal until the end of the second chirp.

Now at the processing part my doubts are the following:

If the two nodes run the same clock i would know perfectly where the signal starts from the end of the first  chirp, after the fixed period of pause. In this case the signal i try to demodulate would have no phase shifting, correct?!

If my clocks run differently, i will reconstruct a different signal, right? If yes, why is shifted in phase? Or what does it mean shifted in phase in such case. Can you point me out how to adjust, if possible, this offset or to know where to sample the signal?  Is this effect something like a Doppler?

Please consider that my demodulation is in passband and i don't know if there exists any methods which would not require baseband transformation.

Thanks in advance

Regards


texts: How to efficiently control an FIR's magnitude response by altering its phase spectrum

texts:
Question:

Extensively searching the space of all possible vectors of length $n$ to satisfy a (non-overdetermined) requirement is possible in principle. Hence, there is a way to calculate a complex FIR kernel of length $n$ whose magnitude response goes through $2n$ given points, when its phase response is allowed to be completely arbitrary. (Corollary, a real FIR kernel with $n$ given magnitude points and $n$ mirrored points can be calculated.) The extensive search method is hopelessly inefficient, though. Is there a practical way to calculate an FIR with such requirements?

Context:

It is common DSP practice to use the DFT to convert an array of real magnitude points $M = \{m_0, m_1, \dots, m_{n-2}, m_{n-1}\}$, $m_i \in \mathbb{R}$ along with an array of real phase points $P = \{p_0, p_1, \dots, p_{n-2}, p_{n-1}\}$, $p_i \in \mathbb{R}$ to an FIR kernel $K = \{k_0, k_1, \dots, k_{n-2}, k_{n-1}\}$. $K$'s elements are generally complex, but if $M$ is symmetric (i.e. $m_i = m_{n-i}$ for $0 < i < \frac{n}{2}$) and $P$ is conjugate symmetric (i.e. $p_i = - p_{n-i}$ for $0 < i < \frac{n}{2}$), then $P$'s elements are completely real.

The magnitude response (absolute value of the DFT) of $K$ will always match $M$ at the individual points $m_i$. Between these fixed points, however, it can attain very unintuitive and undesirable values. This undesirable behavior can easily be visualized by appending a large number of zeros to $K$ and then calculating the absolute value of the DFT. Changing any $p_i$ changes the behavior of the magnitude response of $K$ in between the $m_i$. Thus, two FIRs derived from the same magnitude spectrum $M$, but with different phase spectra $P$, can have very dissimilar magnitude responses. This is illustrated by the figure below:

The upper plot shows three randomly chosen examples of $P$, called $P_0$, $P_1$ and $P_2$. The lower plot shows a randomly chosen example of $M$ and the resulting magnitude responses when an FIR kernel is formed using $M$ as the magnitude spectrum and one of the $P_i$ as the corresponding phase spectrum. $n$ was chosen to be $16$ here. The three magnitude responses all go through the fixed set of points defined by $M$, but vary considerably everywhere else, due to their different phase spectra.

There are many audio applications where we do not really care about the phase response of a filter, but we absolutely do care about the magnitude response. Given a magnitude spectrum $M$, we can hence regard all the $p_i$ as degrees of freedom to shape our magnitude response. More specifically, this should allow us to define $n$ additional magnitude points that the magnitude response of $K$ is required to go through. These additional requirements should be able to help to substantially reduce any undesirable magnitude response behavior of $K$ without the need to increase its length $n$. Calculating $P$ to satisfy such additional requirements seems nontrivial, though.

Are my assumptions correct? Am I missing something? Is there any practical (polynomial complexity) way to do this?

Edit:

This seems to come down to the question of wether a complex polynomial can feasibly be interpolated by a number of absolute value requirements. The coefficients of the FIR kernel $K$ define a polynomial in the complex plane and $K$'s frequency response lies on the unit circle. Any magnitude requirements are positions on the unit circle with a value which the absolute value of that polynomial is supposed to go trough. Consequently, I asked a more specific question on Mathematics Stack Exchange. If I gain any further insight from there, I will update or answer my own question here.


texts: How to Actually Use FFT to Implement Low Pass Filter

texts:
have searched high and low for a practical guide to implement a simple 50Hz low pass filter on a signal.I understand the theory but its the practical but that's missing. If I was using MATLAB it would be trivial because there are so many examples but I am coding this by code.

This is where I am at:

1: I have a series of 8192 bytes sampled at 1000Hz in an array (0-256)

2: I have generated the frequency domain representation of a simple step filter in a symmetric array (+/-F, Ire, Iim) :

0,1,0
0.1220703125,1,0
0.244140625,1,0
0.3662109375,1,0
0.48828125,1,0
0.6103515625,1,0
.....
-0.6103515625,1,0
-0.48828125,1,0
-0.3662109375,1,0
-0.244140625,1,0
-0.1220703125,1,0

3: I have used Inverse FFT of the filter to obtain the time representations of the filter as (t, Ire, Iim) [as expected the im parts are ~0 and the filter gives the expected decaying ringing]:

0.000,819,0
0.001,805.600848087199,4.11892742135933E-12
0.002,766.191083986647,7.80009390410896E-12
0.003,703.07869137741,1.07704956064936E-11
0.004,619.930866652693,1.26768595620774E-11
0.005,521.519198753447,1.32877042702262E-11
0.006,413.390428478234,1.26576527037514E-11
0.007,301.48763797423,1.08400510789863E-11
8,0.008,191.750273447113,7.79765141345479E-12
....
8.186,413.390428478235,7.79385628257856E-09
8.187,521.519198753448,9.84392145575441E-09
8.188,619.930866652693,1.17116023545805E-08
8.189,703.07869137741,1.32913150485692E-08
8.19,766.191083986645,1.44921696865197E-08
8.191,805.600848087196,1.52439926792702E-08

OK so what do I do with this data???

From what I have read in the literature the thoery is that in the freq domain you multiply the spectra but because this is the DFT the usual course is to convolve in the time domain but the question remains:

1: Do I simply multiply the two time series together (doesn't feel right)?

2: If not how is the weighting applied?

3: How do I normalise any output back to 0-256 values?

I realise this is a pretty basic question but there doesn't seem to be a really basic answer or example out there :-(

Thanks

iced


texts: Calculating cross-correlation using Walsh-Hadamard transform

texts:
I am trying to implement MLS method of measuring impulse responses. There is an article describing the method: http://www.commsp.ee.ic.ac.uk/~mrt102/projects/mls.... As I understand, to get an impulse response of a system I need to calculate cross-calculation  of input signal - MLS and output signal of a system.
I tried to do measurements on simple lowpass RC-filter. The result I get is strange, it is a noisy signal that does not look like lowpass filter signal response. I tried to use convolution and circular correlation but it didn't help.
I suspect the problem is related with calculation methods. I calculate correlation using the convolution theorem that allows to calculate convolution using Fourier transforms of signals. But, the article describes Walsh-Hadamard transform to use in the method.
Also I checked the form of MLS, the generated signal is correct, but I am not sure should it be periodic or not.
My questions are:
Can I get wrong results because I don't use Walsh-Hadamard transform during calculations?
Could you, please, give a sample of code in Matlab or Python that calculate correlation using Walsh-Hadamard transform?
Can there be other sources of problems, for example, distortion related with periodicity of MLS-signal?
Thank you in advance.
P.S. I can show screnshots of wrong impulse responses I got if it is required.


texts:
I am trying to make FFT simulation in Matlab by generating noise added
two sinus waves in 60Hz and 100Hz.

After adding the noise into these signals then I have applied the FFT
as I put my Matlab code below.

But I am in difficulty interpreting the FFT plot which shows two spectral peaks in the plot diagram below.

Could you please explain how to interpret this spectrum?
Why we see two peaks in the spectrum? How can we reduce the spectral to show only the frequencies of 60Hz and 100Hz?

%%%% Noise_Added_Two_Sinus%%%
>> f1=60;
>> f2=100;
>> fs=512;
>> t=0:1/fs:2-1/fs;
>> x1=2.4*sin(2*pi*f1*t);
>> x2=0.96*sin(2*pi*f2*t);
>> y=x1+x2+randn(size(t));
>> F=fft(y);
>> plot(abs(F))


texts:
I did a bit of searching on this site but surprisingly, I didn't find much relevant information, and my knowledge of DSP is very very limited.

My goal is quite simple : I want to program an algorithmic reverb in C++ which sounds really good. Or more exactly, the best option would be to let the end user choose the trade-off between quality and cpu-usage.

From what I found out so far is that, to create a reverb, you must feed the dry signal into a early reflections algorithm, then to a late reflections algo. Is this correct ?

Now, I found an extensive article on the late reflections part, using Feedback Delay Network (REDUCING ARTIFICIAL REVERBERATION REQUIREMENTS
USING TIME-VARIANT FEEDBACK DELAY NETWORKS). From what I've read, FDN is a high quality, not too expansive (cpu wise) way of simulating the late reflexions. Furthermore, I guess you can control the quality/cpu charge trade-off by changing the number of delay lines.

However I have absolutely no idea of how to program the early reflections algo (remember ? I'm  really ignorant in the DSP domain).

It would sound logical to me to use a kind of multi delay, which has the advantage to be easy to program and computationally inexpensive. But that sounds too simple to be true.

Furthermore, my instinct tells me that one or several filters must be included somewhere in the signal path.

Could some one please clarify this topic a bit for me ?

Two notes :

I'm not after a convolution reverb at all. I don't really care about the realism of the reverb, but instead I want a good sounding, tweakable, not to cpu-hungry reverb.
Also, the coding part isn't what worries me, I'd ask on stackoverflow otherwise. It's really the DSP part, and just that part, which I'm after :)


texts: Why Successive Approximation, Analog to Digital converter is most commonly used?

texts:
There are various types of ADC in used, like,

ramp-compare ADC
Wilkinson ADC
pipeline ADC
Flash ADCs
Sub-Ranging ADCs
Folding ADCs
Pipelined ADCs
Successive Approximation (Algorithmic) ADCs
Integrating (serial) ADCs

but  Successive Approximation (Algorithmic) ADCs (ADC0808) is most commonly used.
I do not understand the reason as why it is better than any other ADC ?

According to this comparison,

There is absolutely nothing so special about it as compared to other, then why it is most commonly used for education purpose and small electronic projects etc ?


texts: Necessary Conditions for stability in z domain?

texts:
What are the necessary(must) conditions for stability in z domain?
I am sure about one(ROC must include unit circle)
Is there any other such condition which states that there shouldn't be any poles in ROC?


texts: Question about convolution in time and frequency domain

texts:
My question is about the number of data points when doing a convolution (or correlation) in the the time or frequency domain. let's say we have two signals, $a(t)$ and $b(t)$ each of length $16$. When we do the convolution as follows:

$$c(t) =a(t)\star b(t)$$

$c(t)$ consists of $31$ points, right? Now if we transform these into $A(f) = \mathrm{FFT}(a(t))$ and $B(f) = \mathrm{FFT}(b(t))$, $A(f)$ and $B(f)$ each consist of $16$ points still, were the first half is a mirror of the second half. The convolution is just a multiplication in the frequency domain:

$$C(f) = A(f)\cdot B(f)$$

I believe $C(f)$ should be $16$ points. I'm fine with that much, but if we transform $C(f)$ back to the time domain, doesn't  $c(t)$  end up having $16$ points?

I'm sure I'm just not thinking correctly about this. Thanks for any clarification!


texts: Why is scale space (DoG) needed to detect scale invariant features?

texts:
I would have a theoretical question, about the detector part of the SIFT algorithm, which as it is explained by D. Lowe in his paper, used the DoG to detect keypoints in an input image.

Problem formulation: Let $x$ be my original image, and let $y$ be a scaled-by-two version of $x$.

In the SIFT algorithm: I don't understand how the SIFT detector can detect a same given keypoint from image $x$ in image $y$. I don't understand how the "scale space" DoG pyramid is used exactly to achieve this "scale invariance".

Please note that I have watched quite a few videos and read a few tutorials about it. I understand how the DoG pyramid is constructed, by successive blurring of the image. I also understand that keypoints are the maxima or minima from comparing to adjacent images in the DoG pyramid. (I understand that the difference of 2 blurred images with different variances gives an edge-like image.)

However, for me, in every tutorial, they sort of "skip" the fundamental part that would clearly explain why it makes it scale invariant.

In other words, I understand that different gaussian blurs (increasing variance) allow to highlight different types of details in images (blob detector). I'm ok with that, but I don't see any relation or reason intuitively why this has anything to do whith "scale" and why and how that is used to detect a same keypoint in image $x$ and $y$ with for example image $y$ zoomed in twice.

I have essentially the same (unanswered) question as him.


texts:
From acceleration data I made a Power Spectral Density (PSD) and calculated the RMS via:
$$a_{RMS} = \left[\displaystyle\int_{f_1}^{f_2} a_{PSD}(f) df\right] ^{\frac 12}\tag{1}$$
Now, I want to calculate the RMQ (Root Mean Quad) from the PSD. Considering the fourth-power relation for the RMS I think the following formula should translate the PSD to the RMS:
$$a_{RMQ} = \left[\int_{f_1}^{f_2} a_{PSD}(f)^2  df\right] ^{\frac 14}\tag{2}$$
But I did not manage to prove this.
Can someone provide me some extra information/a prove about this problem. I can not find anything about calculating the RMQ from the PSD. Thanks in advance!


texts: Compute Similarity between Fourier Transforms

texts:
I'm looking to compare the Fourier Transforms generated by accelerators and gyroscopes that collected data of people walking. I've looked to see if there is a standard form of comparison, but I have yet to find one.

Here's my specific use case: I'm looking to compare similarity between the dominant signals of two Fourier Transforms (i.e. how close are they to being the same frequency/magnitude). Is there a specific metric that would help identify that?

If not (or perhaps irrespective of my goal), what are some common ways of comparing two Fourier transforms? I've found surprisingly little on the subject.


texts: Pitch detection when successive notes are the same

texts:
I want to program a microcontroller in C but I am asking about the theoretical aspect of it. I think I've got an idea what I am going to do but there is one thing I can't quite get. To begin with, I will use a guitar, so there is a 82Hz-987Hz range. So roughly, I'm going to low-pass the data on 1KHz, "window" them and apply FFT. Supposedly, I am going to have a good estimate.

Since highest frequency is ~1KHz, I decided my sample rate to be 2KHz (haven't tested yet if it works).

Suppose I have a 2 seconds track with 1 note each second. If the first note is C, my result would be 2000 C notes for the first second right? We can say that every time there is a note change, that's when you recognize a note, so if the second note is D, the program will acknowledge C, then D (and not 2000 times C and 2000 times D). The problem I see, is that if the second note is also C, how would I know that indeed there is a note hitting?

I believe I need an extra step to the process but I can't quite find what would that be. Might be that I have to go a different route though.. (if this is the case just mention the method, no need to elaborate on it and lose your time, I will read about it).


texts:
The figure below shows in dashed lines sinusoidal signals of the same frequency at three different phase shifts.  The signals are then sampled such that the sinusoidal frequency is exactly a half of the sampling frequency, i.e. the frequency of all the sinusoids is the Nyquist frequency.  The samples taken from this signal are represented by the circles.

From the figure, it seems that the amplitude of the digital sinusoidal signals is dependent upon the sampling rate and instant of sampling.  In fact if the sampling times coincide with the zero-crossings of the sinusoid, then no signal will be detected at all.

I had initially thought that sampling a bandlimited analog signal at the appropriate sampling frequency would enable perfect reconstruction, but this counter-example has left me stumped.  It seems that this sinusoid will generally not be reconstructed properly if digitized and then reconstructed at this rate.  Have I gone wrong in my understanding, and if so, can someone please point me in the right direction?

Cheers!


texts: What is the correct way to apply multiple IIR filters to a signal in MATLAB

texts:
I am attempting to apply multiple IIR filters (one band-pass, one high-pass, and one low-pass;each with Fs of 44.1 kHz ) to a signal in MATLAB by using the filter function.  When I apply one signal and sound it, I can hear it. But as I start to apply the others, I can no longer hear the signal. I have the following code:

 audio = handles.filename;
 [x,Fs] = audioread(audio);
 lowpass = IIR_Lowpass_Butterworth_CFL1000;
 highpass = IIR_Highpass_Butterworth_CFH10000;
 bandpass = IIR_Bandpass_Butterworth_CFL1000_CFH10000;

 y = filter(lowpass,x);
 z = filter(highpass,y);
 q = filter(bandpass,z);

 sound(q);

Am I applying the filters incorrectly?


texts:
Please, pardon my ignorance. This may be a simple question, but I am still unable to explain what is happening.

I am using a turtlebot to listen to sound emitted from a stationary sound source.

My experiment setup is this:

One stationary sound source emits tone of 500Hz.
A turtlebot robot located about 3.5m away listens and measures the amplitude of the sound source (used FFT and filtered out 500Hz) at a rate of 100Hz while moving linearly toward the sound source at a velocity of 0.1m/s.

The result of the measurements is shown below.
NB: X axis is distance in metres, Y axis is amplitude measured by the turtlebot.

My question is this:
1. Why is there a noticeable dip in the amplitude as the robot moves closer to the sound source?
2. Please can you help me with links to resources that can help me explain this behaviour?


texts: Artifacts in image after super resolution using delaunay triangulation in MATLAB

texts:
i have to do super resolution of two low resolution images to obtain a high resolution image.

1st image is taken as base image and the second image is registered with respect to it . i used SURF algorithm for image registration . A Delaunay triangulation is constructed over the points using a built-in MATLAB delaunay function . The HR grid of size is constructed for a prespecified resolution enhancement factor R Then HR algorithm for interpolating the pixel values on the HR grid is summarized next.

HR Algorithm Steps: 1. Construct the Delaunay triangulation over the set of scattered vertices in the irregularly sampled raster formed from the LR frames.

Estimate the gradient vector at each vertex of the triangulation by calculating the unit normal vector of neighbouring vector using cross product method.Sum of the unit normal vector of each triangle multiplied by its area is divided by summation of area of all neighbouring triangles to get the vertex normal.

Approximate each triangle patch in the triangulation by a continuous and, possibly, a continuously differentiable surface, subject to some smoothness constraint. Bivariate polynomials or splines could be the approximants as explained below.
Set the resolution enhancement factor along the horizontal and vertical directions and then calculate the pixel value at each regularly spaced HR grid point to construct the initial HR image

now i have the results shown below

now for one kind of data set i get this result that has a few pixels black and wite in a random manner for the other type i get thin parallel lines all over image after super resolution the results are attached

for first set random pixels appear on super resolution,for second set random pixels with parallel lines appear, these pixels appear on differnt pattern when i change transformation type from affine to projective in estimate geometric transform command, more over for other data sets there are deformation like parallel lines in edges.for some images it is quite well and does not give any deformation

any one can tell me the reason, i have figured out may be its demosaicing but i am not sure,because i dnt have much understanding of it , moreover can it be a bug in my code but it behaves different for different images, i have increased the size by super resolution twice.


texts: How to detect time of arousal with excitation signal in speech signal?

texts:
I have several seconds of speech signal and I have to detect moments in time of arousal with excitation signal, which is creating by larynx. I'm looking for useful reference since several days. I know that a speech signal is a convolution of excitation signal with a filter response and I've read about cepstral analysis and linear prediction, but I still do not know how to take advantage of these methods for my puspose. I think that I have to separate these two signals in some way.

I found the fundamental frequency by create autocorelation function and reading maxima, which gave me the period of fundamental tone. But actualy I don't know how to use it.
I need to do it in matlab. Please help me.

Here is good example what I need.
http://www.iba.zut.edu.pl/dydaktyka/pliki/referat_js/image002.jpg
In the top plot these impulses are marked as short lines and this is what I need to do.


texts:
I have the Fourier Transform or FFT of a signal, call it $F(f)$. I would like to apply a "frequency-dependent frequency shift" to it. Specifically, frequency $f$ should be shifted to $3f^2$, for all $f$. Physically, this signal is undergoing dispersion.

Say I have a signal of a one-dimensional wave (sound amplitude, water height, or voltage, etc.) as a function of time at a specific location, and I'm trying to estimate the spatial signal at a single instant in time. Let the dispersion relation be $k=3f^2$, where k is the wavenumber. This procedure will change the temporal signal to the spatial signal

This seems like a linear system to me, so how would I find the frequency response $H(f)$ of the system? Knowing the frequency response, I can simply multiply $F(f)$ by the $H(f)$ to find the resulting signal.


texts: Generate binary signal and convert to integer in Matlab

texts:
I am going to test an algorithm and need a test signals.
My first test signal is a binary signal of 1 or -1. I was thinking about how to implement such a sequence in matlab and wrote this:
binary_data = rand(1,M) > 0.5;
b = 2.*binary_data - 1;

My next step is:
dataInMatrix = reshape(b,length(b)/k,k); % reshape data into binary k-arrays, k=log2(M), M -size of signal constelation
dataSymbolsIn = bi2de(dataInMatrix);

error:
Expected input number 1, binary vector b, to be nonnegative.
How to convert such binary sequence to integer if i need create symbols?
EDIT 1.
Guys, I have generated a test signal and it look like this:

i generate a sequence 1 or -1
use a raised cosine filter

How can I get a blue lines...show samples ( digital signal)?


texts: Alternative function for Matlab Stationary Wavelet Transformation

texts:
I'm currently working on a project of mine where I have to use stationary wavelet transformation. I found a perfect way to do the transformation using wavelet toolbox of the Matlab. But since it is not free, I need to find a free function to do the exactly same thing SWT function of wavelet toolbox does.

What I Found

I found a Matlab functions package called WaveLab. They have defined a function called FWT_Stat to do Stationary wavelet transformation.

FWT_Stat -- stationary wavelet transform
Usage
 StatWT = FWT_Stat(x,D,qmf)
Inputs
 x        array of dyadic length n=2^J
 L        degree of coarsest scale
 qmf      orthonormal quadrature mirror filter
Outputs
 StatWT   stationary wavelet transform table
          formally same data structure as packet table
          log_2(n)-D scales by n elements

Below shows information about SWT function from wavelet toolbox.

SWC = swt(X,N,'wname') computes the stationary wavelet decomposition of the signal X at level N, using 'wname'.

Problem

FWT_Stat function asks for 3 parameters and final parameter is orthonormal quadrature mirror filter.But I couldn't find a way to add haar filter as the input for 3rd parameter .If I used SWT function from wavelet toolbox, all I have to do is add wavelet name (e.g: harr) instead of wname and swt function take care of all other things.

So my problem is, how can I add haar filter as the 3rd parameter for FWT_stat function on Matlab without Wavelet toolbox.

Or, Is there any alternative function (which is free) to do the exact same thing as SWT function from wavelet toolbox?


texts: Effect of Gaussian Smoothing to mean image intensity

texts:
I am implementing Perona-Malik Diffusion in using  the formulas:

Creating image at time k+1 using equation

This is my implementation:

function gij=g(duij,lambda)
duij=-(duij.^2/lambda^2);
gij=exp(duij);
end

function deltauij=gradMag(u)
xx=[u(:,2:end),u(:,end)]-[u(:,1),u(:,1:end-1)];
yy=[u(2:end,:);u(end,:)]-[u(1,:);u(1:end-1,:)];
xx=(xx.^2)/4;
yy=(yy.^2)/4;
deltauij = xx+yy;
deltauij = sqrt(deltauij);

function u=nld(u,delt,lambda,sigma)

filt= fspecial('gaussian',[3*ceil(sigma),3*ceil(sigma)],sigma);
utemp=imfilter(u,filt);
deltauij=gradMag(utemp);
gm=g(deltauij,lambda);
uxyplus=[u(:,2:end),u(:,end)];
uxyminus=[u(:,1),u(:,1:end-1)];
uxplusy=[u(2:end,:);u(end,:)];
uxminusy=[u(1,:);u(1:end-1,:)];

gxyplus=[gm(:,2:end),gm(:,end)];
gxyminus=[gm(:,1),gm(:,1:end-1)];
gxplusy=[gm(2:end,:);gm(end,:)];
gxminusy=[gm(1,:);gm(1:end-1,:)];

xx=0.5*(gxplusy+gm).*(uxplusy-u)-0.5*(gxminusy+gm).*(u-uxminusy);
yy=0.5*(gxyplus+gm).*(uxyplus-u)-0.5*(gxyminus+gm).*(u-uxyminus);
ddu=xx+yy;
del = delt*ddu;
u=min(255,max(0,del+u));

And I am calling iteratively  I=nld_2(I,delta_t,lambda,sigma)
I am calculating mean intensity value of image in every iteration.The result is quite weird .Mean vs iteration plot is lik this

Since this is a smoothing operation I am expecting it to decrease mean intensity somehow but what I get is an increasing mean intensity.

Is this normal or I am missing anything?


texts: audio features extraction using restricted boltzmann machine

texts:
I want to extract Audio Features using RBM (Restricted Boltzmann Machine). For this, I am giving the spectrogram (PCA whitened) as an input to the RBM. For each audio file, The spectrogram is a matrix with no. of columns fixed but with different number of rows for each audio file. My question how can I train my RBM, or how can I extract the features from audio using RBM, given this spectrogram matrix. I read in a paper by Honglak Lee, paper title Unsupervised Feature Learning for Audio Classification using convolutional deep belief networks.
http://machinelearning.wustl.edu/mlpapers/paper_files/NIPS2009_1171.pdf
"We then trained 300 first layer bases with a filter length of 6 and a max-pooling ratio of 3."
First, what is meant by bases here. (They have used Convolutional Deep Belief Networks, so I guess, bases do not mean weights here).
Second, what do they mean by using a filter length of 6? How can I do it? Any hint will be appreciated. (I am new to RBM)


texts: I cannot find reference (paper) of this relation $u(t)t^{α} ↔^{FT} f^{-(α+1)}$

texts:
I am struggling a lot in finding the reference (book or, even better, academic paper) that states that:
$u(t)$     $t^{α} ↔^{FT} f^{-(α+1)}$
as I found here An Interesting Fourier Transform - 1/f Noise in the formula in the first box of the link.
I searched a lot, but I did not find any reference, and they are very important to me.


texts:
The measure of a given frequency $\omega$ in a signal $x(t)$ is:

$\frac{1}{N}\sum\limits^N_{t=0}x\left(t\right)e^{^{-i \omega t}}$

This is basically an average of the correlation between the signal and a complex sinusoid rotating at $\omega$ frequency over some time range (N).

My goal is to realize this as a IIR filter (the former definition is virtually a FIR filter), so that it can be implemented as a sort of resonator which is fed an input stream of data, from which the frequency over the last N samples can be measured, continuously?

Currently, I've made this working model:

$y\left(n\right)=e^{-i \omega n}x\left(n\right)$

$X\left(n\right)=y\left(n\right)+e^{\frac{-1}{N}}\left(X\left(n-1\right)-y\left(n\right)\right)$

Where $X(n)$ will be the complex result of the correlation over the last N samples, utilizing a moving average IIR filter. While this works, it is quite inefficient in my current, naive implementation:

for each sample:
    t0 = x * sample;
    real = t0 + pole * (real - t0);

    t0 = -y * sample;
    imag = t0 + pole * (imag - t0);

    t0 = x * c1 - y * c2;
    y = x * c2 + y * c1;
    x = t0;

where x and y corresponds to the cosine and the sine of the complex sinusoid, computed using a IIR oscillator filter. My hunch is that it will be possible to implement some standard filter with an adequate damping such that it will be equivalent to the former equation / code? And basically function like a resonator, where the amount of resonance is equal to the average of the correlation over the last N samples.

I've also been looking at the Goertzel algorithm, and possible ways to use that theory, even though it seems to suffer from scalloping/inability to resolve frequencies not quantized to integers between 0 and N?

Sorry for the redundancy, I'm still very new to filters and DSP in general, I'm just hoping someone can point me in the right direction and tell me I'm not grasping for straws :)

Thanks

Edit - appendix:

goertzel algorithm

Scalar sine, cosine, coeff, q0(0), q1(0), q2(0);

sine = sin(omega);
cosine = cos(omega);
coeff = 2.0 * cosine;

for (int t = 0; t < size; t++)
{
    q0 = coeff * q1 - q2 + data[t];
    q2 = q1;
    q1 = q0;
}

Scalar real = (q1 - q2 * cosine) / (size * 0.5);
Scalar imag = (q2 * sine) / (size * 0.5);

goertzel vs. correlation of sinusoid