texts: Understanding Houghlines command in Matlab

texts:
I created this picture with paint to understand the hough transform but i realized something funny by using houghlines(bw,theta,rho,peaks) in matlab.

I=imread('deneme.png');
I=rgb2gray(I);
bw = edge(I,'sobel');
[H,theta,rho] = hough(bw);
peaks  = houghpeaks(H,1); % I select a peak which is maximum
lines = houghlines(bw,theta,rho,peaks);

Matlab found

peaks=[1097 55]

If we look at the lines variable

And using plot command , i showed where these points provided by lines.

imshow(bw)
hold on
plot(lines.point1(1),lines.point1(2),'r*')
plot(lines.point2(1),lines.point2(2),'r*')

My question is why these red points here.

Instead of being here , should they be at the start of the line and end of the line ?

And where those points come from.

r=x cos(theta) + y sin(theta)

from rho and theta , we can obtain this equation. But not the points

Thank you


texts: How to segment femoral and tibial knee cartilage individualy?

texts: I am researching in medical image processing area. I want to segment knee cartilage into femoral cartilage and tibial cartilage. But two part have same intensity and only different position(femoral cartilage is above and tibial cartilage is below). I am finding the techniques or method to segment them individualy.Could you help me show to me that techniques? Thank you so much![enter image description here][1]


texts:
What is $N$,N-point Discrete fourier transform? Is it different from 2D Fourier transform?
and how to compute $N$, N-point Discrete Fourier Transform of a given laplacian filter kernel? say for a given filter {{0, 2, 0}, {1, -2, 1}, {0, 2, 0}}, where $N = 103$ and frequency and spatial indices ranges from -100 to 100, without using fft2 and ifft2 functions.


texts: multiplying the received signal by carrier, in OFDM system

texts:
Sorry for this question, I'm still new in signal processing area.

In OFDM system, let $x$ be the received signal, what do we call the process of multiplying the received signal $x$ by $e^{-j2\pi F_c t}$ , where $F_c$ is frequency carrier and $t$ is interval 0 : sampling rate: sampling rate * length of signal.

So the question, what is that function? and why do we use it.


texts: Is it possible to "equalise" a signal by deconvolving the impulse response of the room in which it is to be played?

texts:
I am using a test sweep with a flat power spectrum and linear group delay (Optimized Aoshima's Time Stretched Pulse) to measure a room's frequency response. Having obtained the impulse response of the room, I thought I would be able to deconvolve it from my test signal, record the result in the same spot and under the same conditions, and obtain an approximately flat signal. However, the result of the deconvolution sounds like random noise, not at all like my original sweep.

I've tried naive deconvolution by inverse filtering, regularising my IR to avoid zeroes in certain frequencies, even ignoring all phase information in the IR and divinding just the magnitude of the results, all to no avail. In every case, the result of the deconvolution of sweep with the room IR sounds like random noise. Any ideas where I might've gone wrong?


texts:
I want to implement basic noise reduction. I'm recording audio using a microphone mounted on a drone. I would like to filter the audio file from drone noise.
What's the best way to do it ?

My first idea is to first record sound that is emitted by a drone only and convert it using fourier transform to see which frequency has the most energy.
Then I could do the same thing with audio file with drone noise.

My question is how to cancel this drone noise ?

This is the spectrum of the sound emitted by propellers:

Most of its power is in the low frequencies.

Edit.
Community noise, including aircraft and motor vehicles, typically ranges between 50 Hz and 5,000 Hz. So basically if I filter noise from propellers I might also filter sound from vehicles?


texts: DFT vs FFT: does odd number works for DFT only?

texts:
I am always confused between FFT and DFT. In both algorithm, it has been assumed that the signal is periodic (so my understanding is that, if you have a 20 sample point signal, and you do DFT or FFT, the algorithm assumes 20 sample points repeat).

Now here is the question, someone told me that DFT can take odd number of samples, and spectral leakage could be avoid if use DFT.

I disagree with the person because spectral leakage occurs in both FFT and DFT. And I also disagree because DFT and FFT can both take odd number of samples point, except that it would be slower in FFT case. Is this correct?


texts: Computer Vision - Gaussian Pyramid, Scale Space and Pattern Matching

texts:
I am planning to implement a pattern matching algorithm using something like correlation as a matching metric.

I know that the template I am going to use will, if present, have different sizes in the target images. The interval of sizes might be something like :
x 0.5 , x 0.75 , x 1.0 , x 1.25 , x 1.5 , x 2.0 of the original template size.

To speed up the search it is suggested in various articles to build a gaussian pyramid of both the template and the target search image.

However I am not able to relate in any meaningful way the number of pyramid levels to create (for model and for target) , the object scales as described above, and the gaussian filter sigmas (used when creating the pyramid) !

Could someone shed any light on this ?

For example 1 level of the pyramid reduces the resolution of the image by 2 (is this the same as saying that the image has been scaled by 0.5 ) ?

Than what about a pyramid level of -1 (which should make the image bigger x 2.0) : with what gaussian filter should it be treated before upsampling ?

What about object scales (e.g. x 1.25) which fall in-between pyramid levels ?

Thanks to anyone who could provide some insight/references,
    Todor


texts: How does OpenCV find chessboard corners?

texts: What are the internal workings of the OpenCV function findChessboardCorners?


texts: Design causal FIR filter to approximate hilbert transform

texts:
I'm attempting to design an FIR filter that approximates the Hilbert transform in order to get a $90^{\circ}$ phase shift and unity gain. However, I'm having trouble adjusting the filter to make it causal. I've already read this answer, which was very helpful. Here's what I've got so far.
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal.windows import hamming

N = 50
H = hamming(2 * N + 1)

def coeff(n):
    if n == 0:
        return 0
    return H[int(n + N - 1)] * (2 / np.pi) * (np.sin(np.pi * n / 2) ** 2 / n)

def freq(w):
    return np.sum(
        [coeff(n) * np.exp(-1 * 1j * w * n) for n in range(-N, N + 1)]
    )

if __name__ == "__main__":
    freq_vals = np.linspace(0, np.pi, 1000)
    freq_resp = [np.abs(freq(w)) for w in freq_vals]
    ang_resp = [np.angle(freq(w)) for w in freq_vals]
    # plot frequency response
    plt.plot(freq_vals / np.pi, freq_resp)
    # plot phase shift
    plt.plot(freq_vals / np.pi, ang_resp)
    plt.show()

coeff computes an ideal hilbert transform coefficient using the equation
$$
h[n] = \begin{cases}
\frac{2}{\pi}\frac{\sin^2(\pi n/2)}{n} && n\neq 0,\\
0 && n = 0,\\
\end{cases}
$$
which I've taken from Discrete-Time Signal Processing (3e) p.959. I've restricted $n\in[-50,50]$. Additionally, it applies a Hamming window to reduce the Gibbs phenomenon (see the linked answer above).
freq then computes the frequency response for a provided frequency value, using the equation
$$
H(\omega) = \sum_{n=-\infty}^{\infty} h[n]e^{-in\omega}.
$$
Now, when I plot this (plt.plot lines), I get the magnitude and phase response I expect:
mag (x range is 0 to $\omega$, which I've normalized to 1):

phase:

However, I computed the impulse response values for $h[-50],h[-49],\ldots,h[50]$. So, this filter is not causal. I believe this filter should be time-invariant, so I attempted to simply shift it so that $n$ would now range from 0 to $2N+1$. Here's the adjusted freq function:
def freq(w):
    return np.sum(
        [coeff(n - N) * np.exp(-1 * 1j * w * n) for n in range(0, 2 * N + 1)]
    )

This produces the same frequency magnitude response, but with the following "incorrect" phase response:

I expect I've misunderstood something fundamental about this process but I'm not sure what. How can I make this filter causal with the desired phase response? Why does my attempt to make the filter causal keep the correct magnitude response but not phase response?


texts:
I have a question on CP-OFDM.
In recent years people are talking about filter bank approach where they are trying to filter out the out of band side lobes using filters. Besides they are presuming that this approach may provide some robustness and spectral efficiency since they can get rid of CP. Is this true and what aspects we can remove CP and make spectral efficiency equal one.

Ref.
F. Schaich, T. Wild, Y. Chen , „Waveform contenders for 5G – suitability for short packet and low latency transmissions“, accepted for IEEE VTCs’14, Seoul, Korea, April 2014


texts:
We're trying to implement a "simple" speech recognition application in MATLAB (isolated words from a very limited dictionary). We've been trying the following methods:

Extract MFCC coefficients for each frame of the word, and then compare to templates that we have recorded previously using DTW, and take the template that gave the minimal distance as the recognized word. (note: the "templates" are 10 recordings of each of the 6 dictionary-words).
Extract MFCC coefficients for each frame, and run an SVM on the coefficients (for each dictionary word we had a different SVM classifier).
Extract MFCC coefficients for each frame, and then define the feature vector as the vector of all distances from all the templates, then run an SVM on these features.

All three of these methods didn't work well (and gave almost everything wrong results).
We don't know what the problem is.. A few questions that arose, but we don't know the answers, are:

Do we really need 10 samples of each word? or can we "combine" them into a single template? and if so, how?
Should we run the SVM directly on the MFCCs vectors? On the MFCC vector of vectors (for each word)? On the DTW values? or should we combine all the MFCC vectors into a single one?
If we should look for a minimal match for the DTW between the input and the templates, should we give different weights to the templates depending on length, or even account for the differences between the templates (the templates of the word "1" are more distant from one another than the templates of the word "2").

We would appreciate any help, or references to good sources. We thought about using HMMs, but it seems much more difficult and not really possible in our time-frame..

Thanks.


texts: Which one is good Image processing language?

texts: For my academic project, I need to process image to get pixel color and also needs to change the pixel color.I need pattern recognition. It is a web application. All I am trying is to capture the image of debit card by webcam and read the details. Please suggest me programming language.


texts: How do you choose the optimal number of taps in a windowed image scaling filter?

texts:
I'm attempting to use a filter to scale an image and for the most part, I have something that appears to work.  But now I want to be able to generate a filter with an optimal number of taps.  I know in general that the more taps, the better.  But there are diminishing returns in how much each additional tap will increase the quality of the image.  So how would one go about picking a number of taps that makes a good tradeoff between implementation complexity and image quality?  Are there standard measurements that would allow different filters to be evaluated to some "ideal" filter?

Thinking about this more, I think this is really two questions.  One is a matter of deciding which windowed sinc to use (in this case, I've already settled on using lanczos2 or 3) -- that is a more general question that has various tradeoffs depending on the application.  But once that decision is made, the number of taps is easily determined.  For lanczos2, the ideal number of taps is 4 and for lanczos3, the ideal number of taps is 6.  The reason it's not 5 or 7 is that one of those taps will always be zero due to the windowing.

In the more general case, I think the easy way to state how many taps to have for a given window is to simply say:

taps = max - min

In the case of lanczos2, the max is 2 and the min is -2.  Therefore, taps = 4.


texts: How to determine precision needed for sin approximation used for sound synthesis?

texts:
I am doing some procedural sound synthesis in Java. I want to have a sine wave as one of the possible basic sounds. When experimenting with that, I have found the default Java Math.sin implementation causes significant performance issues to me. I would like to replace it with either a table lookup or a simple polynomial (Taylor or Chebyshev) series.
My question is: how to decide what table size or what polynomial degree to use so that the resulting artifacts are inaudible?
I can (and will) perform some basic experiments, but I do not trust my ears and speakers that much.
The sound produced will be a tone with a frequency in range 100 Hz - 10 kHz. I want the result to be indistinguishable from a real sine wave by a human. My sound output will be 16 bit at 44.1 kHz or 48 kHz, if that is important.

how many entries in the lookup table do I need, assuming the table will be used for -pi/2 to pi/2 range?

If I want to use a polynomial series (Taylor or Chebyshev), what degree does the polynomial have to be?


texts: Alternatives to offline Kalman filtering

texts:
Recently I got into vehicle models and filtering in general and immediately faced with the following question.
I have the recorded GPS data from car driving on a highway. However, there is a substantial level of noise present, so I decided to filter this $(x, y)$ data. Currently I am using Kalman filtering (the extended one - bicycle model of the vehicle is assumed) for this purpose, but I somehow have a feeling that Kalman is more suited for online estimation/filtering. I have also tried various other approaches, including median and Savitzky–Golay filters, but they do not take dynamic model into the account which I think is the important information to use while filtering this data.
So, my question is: Is there an alternative to Kalman filter, where model dynamics can be specified, which is more suited for this purpose?
I just once again want to emphasize that I have solely measured time sequence (fixed interval) with $x, y$ coordinates and data needs to be filtered after drive is over (offline usage).
Thank you for your help and I apologize if my question is not stated well - this is my first post:)


texts: Advantages and disadvantages of multiplying sparse vector with sparse matrix

texts: Assume I have a sparse vector $x$ whose few values are non-zeros, that vector is multiplied with sparse unitary matrix $Y$ giving $z = Yx$. Are there any advantages of detecting the sparse vector $x$  based on $z$ the in that case?


texts: Significance of steps involved in flipping process of an image?

texts:
The above image on the right was obtained by following steps

Multiplying the image on the left by $(-1) ^ {(x+y)}$
Computing the DFT
Taking the complex conjugate of the transform
Computing the inverse DFT
Multiplying the image on the left by $(-1) ^{ (x+y)}$

What is the significance of each of the above 5 steps? Here by significance, I mean: what is the benefit of that step? And what loss can occur if we skip that step?

Please note that Steps 1 and 5  involve a multiplication by $-1$ and that $-1$ is raised to the power $(x+y)$.


texts: Zero Padding in Frequency Domain (Spectrum): How to Split the Middle Bin (Nyquist Frequency)?

texts:
When a complex digital signal is converted to its spectrum through an FFT, the result will contain a series op positive and negative frequencies. That is, when the spectrum has N bins. Bin 0 will be the DC, bin 1->N/2 will be the positive frequencies and bin N-1 -> N-N/2 will be the negative frequencies. In this case it is clear that bin N/2 represents the radial frequencies $\pi$ and $-\pi$ at the same time.

Now, If we zeropad the middle of the spectrum (e.g: when upsampling) such that its new size will be M, we obtain a new one as follows. $X$ is the small spectrum. $Y$ is the zero padded larger spectrum.

$X_{[0:N/2[}$ is mapped to $Y_{[0:N/2[}$
$X_{]N-N/2:N[}$ is mapped to $Y_{]M-N/2:M[}$

In the above, I skipped the middle bin ($X_{N/2}$) because it is not yet clear whether it should assigned to the positive or negative section of the new signal.

To better understand what would be the best result, I performed a small test by creating a random spectrum of length M, then converting it to the timedomain, picking out every M/N-th sample and then converting it back to a spectrum, but this time of length N. The result showed that the following is true for any random spectrum:

$X_{N/2}=Y_{N/2}+Y_{M-N/2}$

Further testing revealed that indeed, we can split the value of the middlebin however we want, as long as their sums match the original value. Yet, the interpolations between the samples (when converting $Y$ back to the timedomain) oscillates differently for each particular choice.

The following three images show the effect of the middle bin on the phase of the upsampled signal. The first one places the middle-bin at N/2 in the new spectrum. The second image splits the content over bin N/2 and M-N/2 and the last one maps the middle bin solely to M-N/2

As you see the interpolation goes each time exactly through each of the blue values. The only difference is how the interpolation oscilates.

Is there any rule-of-thumb (and why), or any argument/standard practice on how to split the value of the middle bin when zero padding ?


texts: How to represent a given equation more clear, profressional and short form?

texts:
I have a image $U_{m \times n}:\Omega \to \mathbb R^2$, the output $P$ can be define as
 $$P=\mu J_{m \times n} - U$$
where $\mu = \max \{ u_{ij} : 1 \leq i \leq m, 1 \leq j \leq n\}$, $J_{m \times n}$ be the $m \times n$ matrix whose $i, j$th component is $1$: that is, the all-ones matrix. (This notation isn't quite standard, but it's as close to standard as I know. $J$ is often the all-ones matrix)

However, it is so many sentences for expression the above equation. Do we have more short and standard way to represent it? As my found, the $J_{m \times n}$ can be expressed by indicator function such as

$$P=\max (U)\times 1_{\Omega}-U$$

where $1_{\Omega}$ is Indicator function

Does it equivalent with original meaning? If not, please give me a standard and common expression in image processing. Thanks


texts: Amplitude Modulation with large carrier and suppressed carrier

texts:
From what I have read so far Amplitude modulation with large carrier is favaroble because envelope detector can be used at the receiver to derive the message signal and therefore another mixer is not needed at the receiver.

According to my text book one way to find the envelope of a signal is to lowpass filter the absolute value of the modulated signal.

Using low pass filter on the frequency spectrum $\mathcal F\left\{\lvert (m(t)+1)\cos(\omega t)\rvert\right\}$ will give $m(t)+1$ which is the envelope offset by constant unity.

My question is why can't we use the same technique for AM with suppressed carrier modulation?
Wouldn't low pass filtering of the spectrum of $\mathcal F\left\{\lvert m(t)\cos(\omega t)\rvert\right\}$ give us $m(t)$?


texts: Quadrature component - condition on equivalent lowpass signal

texts:
I have a question regarding the equivalent lowpass signal $x_{LP}$ of a bandpass signal $x(t)$ and the quadrature component.

I'd like to find a condition on $x_{LP}(t)$ and $X_{LP}(f)$ such that the quadrature component of $x(t)$ is zero. (for example $x(t)=a(t)cos(2 \pi f_c t)$).

I know that
$x(t)=Re\{(x(t)+j\hat{x}(t))e^{jw_0t}\}\\
= x(t)cos(w_0t)-\hat{x}(t)sin(w_0t)$

So my idea was that the hilbert transform of x(t) should be zero, in order that the signal has no quadrature component. That would mean that $X_{LP}(f)$ has even symmetry?

I hope that someone can help me with this question!


texts: Analysing mp3 quality with spectrum, how to interpret the results?

texts:
When comparing spectrum images between two different audio sources in order to determine which has better quality, I sometimes come to an impasse. For example, consider the following images from two 320kbps mp3 files:

Optical difference of RGB values:

(here is the difference between waveform-synced versions of the above examples, but it's quite the same)

The second version has a cutoff at 20 KHz, while the first contains some audio information above that frequency, albeit very low in intensity, and probably inaudible.

On the other hand, this extra audio portion at the top of the plot appears to be noise, since it's composed by an almost continuous "blurry" zone. In addition, looking carefully at the two images, some subtle differences can be observed at the lower frequencies.

I'm inclined to think that the second version is better quality, theoretically, since the lower cutoff limit may allow for the encoder to make better use of the available bandwidth, so that it can more acuratelly represent the lower frequencies. That aside, I really don't know how to interpreted those subtle differences between the two images in order to make an objective decision.

So, when comparing spectrum images, what specific details can help determine which encoder did a better job?


texts: Prove that time-reverse yields inverse filter for deconvolution?

texts:
In reading literature on the construction of impulse responses from sine sweeps (e.g. papers by Farina), I see it stated over & over that the way the impulse response is constructed from a (exponential) sine sweep is to take the original sine sweep signal $S(t)$, time reverse it ( $S(T-t)$ ), scale it somewhat in magnitude,  and convolve this with the recorded signal to obtain the impulse response of the (time-invariant, linear) system such as a room's reverberation..  i.e., the claim is that $S^{-1}(t)$ is proportional to $S(T-t)$.

Can anyone explain how to prove this, or point me to a reference that actually shows a 'proof' of this?  The closest I've seen is "it can be shown that...", and a few hand-waving comments about how time-reversal inverts the phase -- but it's not obvious that this 'does the job' for producing a proper inverse from a function with nonlinear time dependence.

I've been trying to show this explicitly (in the continuous-time case) for about 7 hours now, and have a few pages of integrals written out, even with a little help from WolframAlpha :-(.   The exponential term in the frequency makes it hard to compute the Fourier Transform explicitly (you get Gamma functions with two arguments), let alone try to take the reciprocal of the FT and then take the Inverse FT of that.

Thanks!

This is not a homework problem, I just want to understand why it works.  (I know it works because I've coded it up, but it still seems kind of like magic.)

PS- Thought I'd try a linear sine sweep as it might be easier, but you also get a frequency-squared term which again makes it...not at all obvious...that taking $t$ -> $-t$ is going to invert all the phase terms.

P.P.S.- Here are my notes so far; the calculation starts on page 4.  Dropbox link to PDF


texts: Pole Zero plot given a Transfer function

texts:
I've been looking at how to plot zeros/poles based on a transfer function. I found a couple of Tutorials online.

In the first youtube tutorial, the author brilliantly explains how to plot the zeros/poles.

In the second tutorial, the author explains how to obtain the magnitude characteristics of the frequency response, which is also fantastic. In the second tutorial, under the Example in "3.2.3 Transfer function of discrete-time systems", we have the following transfer function

Based on the transfer function, the poles and zeros can be defined as,

a = [1 -2.2343 1.8758 -0.5713]
b = [0.0088 0.0263 0.0263 0.0088]

This is where my confusion starts. based on the first tutorial, i'll have to plot all the zeros/poles along the x-axis (Or am I mistaken?).

But based on the MATLAB command to plot pole and zeros,

zplane(a,b)

I get this plot

The plot for the poles and zeros are scattered all over. How can i plot the poles and zeros manually in the Z-plane given my poles and zeros and obtain a similar output to MATLAB ?

Thanks for your help.


texts: Linearity and time-shifting of $\mathcal{F}\{0.8^n\cos(0.1πn)u[n]\}$

texts:
To preface, this is not a homework related question but purely for self-study purposes.

Hi there, I try to calculate $\mathcal{F}\{0.8^n\cos(0.1πn)u[n]\}$ by using the properties of Discrete time Fourier transform.

Based on the linearity property, we have:

\begin{align}
\mathcal{F}\{0.8^nu[n]\cos(0.1πn)\}&=\mathcal{F}\{0.8^nu[n]\}\cdot\mathcal{F}\{\cos(0.1πn)\}
\end{align}

I got

\begin{align}
\frac{1} {1-0.8e^{-jw}}\pi( \delta(w+0.1\pi)+\delta(w-0.1\pi))\tag1
\end{align}

But if I applied the time-shifting property:

\begin{align}
\mathcal{F}\{0.8^n\cos(0.1πn)u[n]\}
&= \mathcal{F}\left\{0.8^n\frac{e^{j0.1\pi n}+e^{j0.1\pi n}}{2} u[n]\right\}\\
&= 0.5\left(\frac{1}{1-0.8e^{-j(w+0.1\pi) }}+\frac{1}{1-0.8e^{-j(w-0.1\pi) }}\right) \tag2
\end{align}
I am not quite sure if equation (1) and (2) are equivalent to each other. It looks like we have $\frac{\pi} {1-0.8e^{-jw}}$    at (1) when $w$= $\pm0.1\pi$, but I cannot get the same value at (2).

Why I get the different result here? Am I doing something wrong?


texts:
1.If frequently response of a LTI system  exist and  finite then  can we say that, for a periodic or non periodic input signal,  output is also finite?
2.if responses(output) of LTI
system for sinusoids of all the frequencies(inputs),  exist and finite at steady state then we can say that it's frequency response must exist whether system is stable or unstable?


texts: Architecture for high sample-rate FIR/IIR filter on FPGA or ASIC?

texts:
Assume we have a hardware (FPGA or ASIC) which can operate at a fixed processing frequency of $100\textrm{ MHz}$. Assume we have an input data stream of $1\textrm{GSPS}$ (which has been sampled using any SERDES interface, but here the sampling technique does not matter). We would like to implement a filter which can NOT be implemented by a poly-phase decimation filter; for example assume we want to implement a band-pass filter which passes the frequency range $200-300\textrm{ MHz}$ in our example.

What architectures can be used for this purpose?
Can anyone provide a reference for designing such filters?


texts: Improve spectrogram appearance using windowing functions MATLAB

texts:
I have a program which plots data obtained in a NI USB DAQ realtime to a spectrogram plot using the surf function. I'd like to improve the appearance of the plot using some sort of windowing function. My questions are: Which windowing function is the best to use for this? and How do I implement it?
Here's a sample of my current code:

fprintf('Initializing DAQ\n');
s = daq.createSession ('ni');
addAnalogInputChannel(s, 'Dev1','ai1', 'Voltage');

s.NumberOfScans = Fs*Record_Time;
s.Rate = Fs;
s.Channels.Range = [-1 1];
s.DurationInSeconds = 5;

Spect_nFFT = 2^8;
num_overlap = 19*round(Spect_nFFT/20);
start_sample = 1;

DataRdy = true;
% Run in background session

Lh = s.addlistener('DataAvailable', @(src,event) CaptureData(event.TimeStamps, event.Data));
s.NotifyWhenDataAvailableExceeds = 100000*.1;
s.IsContinuous = true;
s.startBackground();
NFFT = 2^nextpow2(L);
P = zeros(8193,50);
F = Fs/2*linspace(0,1,NFFT/2+1);
g=0;

while g==0

    while DataRdy == 0
        pause(.01)
    end
    pause(.5)

    Min_Threshold = 10^12;

    start_sample = 1;
    NFFT = 2^nextpow2(L); % Next power of 2 from length of myRecording
    FFT_Out = fft(myRecording(start_sample:end),NFFT)/size(myRecording(start_sample:end),1);
    FFT_Out = FFT_Out(1:NFFT/2+1);

    A = sqrt(FFT_Out.*conj(FFT_Out));
    P = [P(:,2:end), sqrt(FFT_Out.*conj(FFT_Out))];
    P(P<=max(max(P))/Min_Threshold)=max(max(P))/Min_Threshold;
    surf(1:size(P,2),F,10*log10(P),'edgecolor','none'); axis tight;

    % Create updating spectrogram
    view(0,90);

    xlabel('Time');
    ylabel('Hz');

    drawnow;

    DataRdy = 0;
    uicontrol('Style', 'pushbutton', 'String', 'Exit',...
    'Position', [20 20 50 20],...
    'Callback', 'g=g+1');

end
delete(Lh);
stop(s);

Thanks for any help!


texts: Using Doppler velocimetry with digital signals such as Bluetooth

texts: I have found several sources describing how to use Doppler velocimetry/anemometry to measure the Doppler shift with analog signals, but i cannot find much about using this method for digital signals such as Bluetooth.  I was wondering if it were possible to use Doppler velocimetry/anemometry with digital signals? If so i was wondering if I could use it with the Bluetooth capabilities(2400 MHz – 2483.5 MHz) of an iPhone?


texts: Lower digitisation rate for long distance communcation

texts:
Why is it that usually lower digitisation rates are used for long distance communications? e.g. sound is highly distorted in undersea communications, is this because of the lower digitisation rate used?

why lower digitisation rate is preferable?


texts: Line based Extended Kalman Filter for localization : Measurement Jacobian

texts:
I have to estimate location w.r.t a given map of the environment. The state is expressed as

$ X = \begin{bmatrix} x \\ y \\ \theta \end{bmatrix}$

Given ,

$ X_{t-1} = \begin{bmatrix} x_{t-1} \\ y_{t-1} \\ \theta_{t-1} \end{bmatrix}$ and it's covariance $P_{t-1}, $ a 3X3 matrix .

I have a measurement model which detects lines representing obstacles , so , in one measurement there can be any number of lines , these lines are represented as $$ Z = \begin{bmatrix} \alpha \\ r \end{bmatrix}$$ normal form of line.

I have successfully completed all the steps until the calculation of innovation covariance i.e.

$ S = H_t\hat{P_t}H_t^T + R_t$ ,

where $H_t$ is the measurement Jacobian , $\hat{P_t}$ is the estimated motion covariance and $R_t$ is the measurement noise matrix.

When I have single line measurement , the measurement Jacobian, $H_t$ is a 2-Dimensional $2\times 3$ matrix and the measurement noise matrix , $R_t$ is $2\times 2$. But , when there are multiple line measurements , for example 8 lines , $H_t$ becomes $2\times 3 \times 8$ three-dimensional matrix and $R_t$ a $2\times 2 \times 8$ three-dimensional matrix.

Because of this change in dimension I cannot carry out the calculation of innovation covariance $S$ as $\hat{P_t}$ is a $3 \times 3$ matrix.

How can I incorporate multiple line readings and suppress the three-dimensional matrix to match innovation covariance equation ?


texts:
I'm searching for a method that can smooth a 3D volume whilst preserving the edges in my volume.
I researched anisotropic diffusion filtering and bilateral filtering, but I'm having trouble to evaluate if they can be efficiently used on +30M voxels volumes.
So I need:

Relatively fast algorithm even with 512512200 volumes
Good edge preserving properties
Any help to understand those methods.

Thanks


texts: Algorithm to effectively find perspective transform in motion estimation

texts:
I've shot a series of photographs in a sort of stop-motion fashion, and although the camera was mounted on a tripod, there are tiny motions between each two successive photos.

After some experimentation, I figured that using a perspective transform would be appropriate to apply a correction to a photo so that the difference frame is minimised. This transform basically just moves the four corners of a picture by a given offset, allowing for a trapezoid distortion.

To do this process by hand is impossible in this case, as it takes me like half an hour for each frame. So I am looking for an effective algorithm to estimate the transform parameters (the four corner offsets). I would give the algorithm the maximum magnitude of the offset, and it should then find, as quickly as possible, the parameters (eight values).

Is such an efficient algorithm known? The photos are in colour but could be converted to grey scale, although the lighting changes also over time, so taking a difference image and summing the pixels is not enough, but probably one needs to auto-adjust the contrast and brightness as well.


texts:
The question: Consider transmitting the results of $1000$ flips of an unfair coin where the probability of heads is given by $p_H$. The information contained in an unfair coin flip can be computed:
$p_H\log_{2}(1/p_H)+(1−p_H)\log_{2}(1/(1−p_H))$
For $pH=0.999$, this entropy evaluates to $.0114.$ Can you think of a way to encode $1000$ unfair coin flips using, on average, just $11.4$ bits? (question from https://web.mit.edu/6.02/www/f2011/handouts/2.pdf)
My wrong answer: I thought that I could encode the location of the bits which turn up tails. Since there are 1000 flips, I could encode every flip using 10 bits ($2^{10}=1024$). taking the average expected length to encode each flip and then multiplying by $1000$ for all the flips gives:
$1000[(0.999)(0)+(0.001)(10)]\\
1000(0.001)(10)\\
10$
But I know that any encoding which averages a smaller length in bits than the entropy must have some ambiguity in the message, so since $10<11.4$, what information is my coding system missing?


texts:
Refrence wikipedia article for upsampling and wikipedia article for downsampling

Upsampling , interpolation and expansion are described as synonyms

While

Downsampling , decimation and compression are described as synonyms

Are they mentioned properly as synonyms? Are those 3 terms referring to same operation?


texts: How to plot graph of $e^{-t}$ in frequency domain. What would be the axis? If its Fourier transform is $1 /(1+j\omega)$, then how can we plot imaginary on frequency domain (amplitude vs frequency graph) .


texts: How to remove/reduce noise by using the Morlet wavelet?

texts:
I'm trying to use Morlet wavelet to reduce noise in my signal.
I found formula for Morlet wavelet

$$\Psi(x) = \frac{1}{\sqrt{\pi\cdot \textrm{bandwidth}}} \cdot \exp\left(i \cdot 2\pi \cdot \textrm{centerFreq} \cdot x\right)\cdot\exp\left(\frac{x^2}{\textrm{bandwidth}}\right)$$

(sb told me I shoud use formula for complex Morlet wavelet)

I also got CWT formula

$$C = \int f(t) \cdot \frac{1}{\sqrt{\textrm{scale}}} \cdot \Psi\left(\frac{t-\textrm{shift}}{\textrm{scale}}\right)dt$$

I convert both formulas to java code

public double morletRe(double x)
{
return (1/Math.sqrt(Math.PI*fb))*Math.cos(2*Math.PI*fc*x)* Math.exp(x*x/fb);
}

public double[] cwt(double[] data,double scale, double position)
{
double[] newData = new double[data.length];
    for (int t = 0; t<data.length; t++ )
    {
    newData[t] = data[t]*(morletRe(t-position)/scale)[0])/Math.sqrt(scale);
    }
return newData;
}

But in the results I get a strange chart which is completely not related to my data,
is anyone can point me what I'm doing wrong?
regards.


texts: Envelope of a signal via Convex hull / Edge detection

texts:
Is there a method in order to compute some sort of edge detection / convex hull of a signal x[n] (in orange here). Here of course the orange signal is not convex nor concave, but only piecewise...

(I would like to know if this could be an alternative way of computing the envelope of a signal, alternative to the traditional method input -> abs() -> low pass filter)


texts: gaussian kernel and bionomial coefficients

texts:
The following image is from robert collins' ppt from http://www.cse.psu.edu/~rcollins/CSE486/

How to calculate sigma ?

It shows [1 4 6 4 1] is corresponding to sigma=1.

I learned that usually 3*sigma is best to represent. So kernel size is 7*7 is suited. But its size is 5*5.

I don't know why [1 4 6 4 1] is matched to sigma=1.


texts: what is the theory behind Butter function?

texts:
I want to understand why this code used tan and cos ?

function [ b, a ] = butterTwoBp( dt, fl, fu )
q=pi*dt*(fu-fl);
r=pi*dt*(fu+fl);
N = (tan(q)^2) + sqrt(2)*tan(q) + 1;
M = (tan(q)^2) / N; %M after N because it depends on N
O = -cos(r) * (2*sqrt(2)*tan(q) + 4) / ((cos(q))*N);
P = (-2*(tan(q)^2) + ((  (2*cos(r))   /  (cos(q))   )^2) + 2 )  /   N;
Q = cos(r)*(2*sqrt(2)*tan(q) - 4)/(cos(q)*N);
R = (   (tan(q)^2) - sqrt(2)*tan(q) + 1   )  /  N;

b=[M 0 -2*M 0 M];
a=[1 O P Q R];

I got this code from https://stackoverflow.com/questions/17340198/what-is-the-command-for-butterworth-bandpass-filter?answertab=active#tab-top


texts: The relationship between zero level set and thickness information

texts: I have one problem that need you help. I am try to find the relationship between active contour and thickness information. For defination of thickness we can find at thickness measuring and active contour at chan-vese model. Could you sugess me method or idea to find the thickness of active contour (is called phi) as show in figure


texts: Phase invariant detection of low frequency spikes

texts:
I am a physics graduate student studying a new physical device via simulations. Thanks to SE DSP, reading Lyon's Understanding DSP as suggested (not all the way through it yet though), and a bunch of work, I have been learning, and enjoying, DSP. I’m at a point where I am completely stumped by something, and need some help.

For this project, signal processing appears very simple, as seen in this diagram:

However, what happens inside the "Physics Device" somewhat complicated. Simply put, when the external signal $f(t)$ and the physics device interact, there is some destructive interference. To me, this means there will be a narrow low frequency spike at the time where there is interference. When there is no destructive interference, there is a lot of stuff going on. The outputs relating to the physical device are very ugly compared to a communication system; there are many harmonics and other things happening.

When the physical device I'm simulating is running, and an external microwave signal $f(t)=A\cos(\omega_1 t +\phi_1)+B\cos(\omega_2 t + \phi_1)$ is introduced, destructive interference creates some low frequency oscillations at particular locations:

Input #1: $f(t)=A\cos(\omega_1 t +\phi_1)+B\cos(\omega_2 t + \phi_1)$

Output #1:

Th upper graph shows raw data, and the lower shows this after I apply a Low pass filter. If I use a low-pass FIR filter near DC, it gives me nice peaks to identify where the low frequency oscillations occur.

In these, the first peak near $60\textrm{ ns}$ is related to the first term, and the second smaller peak near $100\textrm{ ns}$ is related to the second term. Also, $B=0.2A$, so the second peak should be shorter than the first peak.

If instead of phase shift $\phi_1$ in the first term, the signal has a phase shift of $\phi_2$, or as: $f(t)=A\cos(\omega t + \phi_2) +B\cos(\omega_2 t + \phi_1)$, we get similar looking raw data:

Input #2: $f(t)=A\cos(\omega_1 t +\phi_2)+B\cos(\omega_2 t + \phi_1)$

Output #2:

But after the low pass filter, the peaks are not-so-beautiful.

Is there an algorithm that will allow me to determine when an external signal is introduced into the physical system regardless of the initial phase shift of the external signal? My eyes can identify from the raw data where the peaks should be, but a low pass filter doesn’t do the job consistently.

I have tried using signal energy, as $\int_t^{t+T} \lvert x(t’)\rvert^2 dt’$, but this fails in two regards.

It does not see the second peak
And it requires some threshold, which doesn’t work well as $A$ and $B$ can differ by several orders of magnitude.

Is there a good algorithm to get nice thin peaks from both signals? Even more, is it something that can be built with existing technology?

I can post data for the above graphs. How can I post it? They are small lists, with ~ 100,000 elements.


texts: Help detecting target in pulsed radar signal

texts:
This signal is from the Pasternack $2.4\textrm{ GHz}$ development kit.  The saved signal appears to be the pulsed waveform.  It is sampled at $20\textrm{ kHz}$ for $4130$ points.  I recorded a background signal, then aimed at a moving target ($\approx 20\textrm{ mph}$ vehicle).

My goal is to remove the background clutter and determine the velocity and range of the vehicle from this data.  So far I have had very little luck.  Here is what I have done:

Get the FFT of background and target signal
Remove the background clutter by subtracting in the frequency domain
Inverse the FFT
Split the signal into $n$ windows
Do the FFT of each and compare to find differences between each window

This method gives me velocities in the $100\textrm{ mph}$ range.  I also tried taking the difference between the background max frequency and target max frequency and still have incorrect velocities.  Are there any other methods to analyze this data?  I feel like the data as provided is useless because of the downsampling and also no time/delay info is included.

A copy of the raw data can be found here.  I think the PRF is about $160\textrm{ Hz}$.


texts:
This post is basically a follow-up to my recent question, which I asked here How to deal with phase incoherence issues after manipulating a signal in the frequency domain?

In short, I'm modifying the STFT and want to resynthesize a waveform that most closely resembles the modified STFT.

After some research, I found out that the standard way to do this is to implement the Griffin Lim algorithm (or one of its variants). The original paper
can be found here https://pdfs.semanticscholar.org/14bc/876fae55faf5669beb01667a4f3bd324a4f1.pdf.

While the algorithm doesn't seem too complicated, I don't think I understand it fully and consequently can't implement it.

I'd appreciate it if someone can have a look at the article and try to explain equation #6 (which is what I need). Some pseudocode/example would be great.

Thank you


texts: What's the expression of the noise for a RF signal?

texts:
I am confused about how to give the expression of noise. Can someone tell me if my equations are correct?

If the received RF signal is expressed as:

where

If I want to add a AWGN, can I give the expression of noise as the following?

If not, what is the correct equation of the noise? Also, in the MATLAB simulation, I often add noise to a baseband signal by using: sqrt(NoisePower) * randn(). What should I do in the case of bandpass signal? Thanks guys!


texts: Clock recovery using Mueller and Muller adds noise affecting EVM or SNR (Two cases - GNU Radio & Python code)

texts:
The implementation of the Mueller and Muller clock recovery that I am seeing adds a large amount of 'noise' or degradation to the Error Vector Magnitude (EVM) to the symbols on the constellation. This is done independently in GNU Radio and a manual script in Python - giving similar results. Is this normal?
I understand that I need to adjust the parameters in this algorithm but I cant get anywhere near an optimisation that is adequate, its too noisy.
In both implementations GNU radio and python , the only impairment to the original data signal is the performance of the RRC pulse shaping filters. You can see that's its an almost perfect received constellation with an EVM < -60 dB using a filter taps of 300 and roll off factor of 0.35. See below image.

However, once I add the M&M algorithm I have no luck getting anywhere close to this and it generates a lot of noise. See below image.

It somewhat makes it difficult to assess the true extent of impairments I may later add in the model, like AWGN. Because if I add AWGN to this then the EVM or SNR just gets worse when the M&M is being used.
Here is the GNU radio flow.
EDIT Note that in my python code the Rx RRC filter is performed before the M&M, but in GNU radio I do it before an RRC because having an RRC before the MM block doesn't work give any result in the flow below.

Below is the M&M python algorithm, its from pySDR.com
samples = input_signal
interp = 16
samples_interpolated = signal.resample_poly(samples,interp,1)
mu = 2 #2 # initial estimate of phase of sample
out = np.zeros(len(samples) + 10, dtype=np.complex)
out_rail = np.zeros(len(samples) + 10, dtype=np.complex) # stores values, each iteration we need the previous 2 values plus current value
i_in = 0 # input samples index
i_out = 2 # output index (let first two outputs be 0)
while i_out < len(samples) and i_in < len(samples):
    out[i_out] = samples_interpolated[i_in*interp + int(mu*interp)] # grab what we think is the "best" sample
    out_rail[i_out] = int(np.real(out[i_out]) > 0) + 1j*int(np.imag(out[i_out]) > 0)
    z = (out_rail[i_out] - out_rail[i_out-2]) * np.conj(out[i_out-1])
    zz = (out[i_out] - out[i_out-2]) * np.conj(out_rail[i_out-1])
    mm_val = np.real(zz - z)
    mu += sps + 0.3*mm_val
    i_in += int(np.floor(mu)) # round down to nearest int since we are using it as an index
    mu = mu - np.floor(mu) # remove the integer part of mu
    i_out += 1 # increment output index
out = out[2:i_out-1] # remove the first two, and anything after i_out (that was never filled out)


texts: Update: Fourier Transform of a shifted and scaled $\operatorname{sinc}$ signal

texts:
Let $x_N$ be the function given by
$$x_N(t)=A\frac{\sin(M\pi(t-N))}{\pi(t-N)}$$
The Fourier Transform of $x_N$ is
$$\begin{align}
X_N(j\omega)&=\mathscr{F}\{x_N\}(j\omega)\\\\
&=\int_{-\infty}^\infty x_N(t)e^{-j\omega t}\,dt\\\\
&=Ae^{-jN\omega }\int_{-\infty}^\infty \frac{\sin(M\pi t)}{\pi t}e^{-j\omega t}\,dt\tag1
\end{align}$$
Enforcing the substitution $t \rightarrow t/M\pi$ reveals
$$\begin{align}
X_N(j\omega )&=\frac{Ae^{-jN\omega }}\pi\int_{-\infty}^\infty \frac{\sin(t)}{t}e^{-j(\omega/M\pi) t}\,dt\\\\
&=\frac{Ae^{-jN\omega }}2\left(\text{sgn}(-\omega/M\pi +1)-\text{sgn}(-\omega/M\pi -1)\right)\\\\
&=\begin{cases}Ae^{-jN\omega}&,|\omega|<M\pi\\\\0&,\text{elsewhere}\end{cases}
\end{align}$$
where we used the Fourier Transform of the sinc function, $\operatorname{sinc}(t)=\frac{\sin(t)}{t}$
$$\mathscr{\operatorname{rect}}(\omega)=\begin{cases}\pi&,|\omega|<1\\\\0&,\text{elsewhere}\end{cases}$$

Let $\displaystyle x(t):=\frac{\sin(2\pi(t-1))}{\pi(t-1)}$. Then we have $A=1,M=2, N=1$ so we get :
$$
X(j\omega)=\begin{cases} e^{-j\omega}&\text{if $|\omega|<2\pi$}\\
0&\text{if otherwise}\end{cases}
$$

Is this example correct? because in my book it says the answer is $e^{2\omega}$

I would hope for someone to assist me in determining if this is correct or not. I have spent entire day on this problem. I would much appreciate any help and thank you :)


texts: I am learning Fourier Transform from many days but till now I am not able to understand what does phase angle image show us or tell us? They say that MAGNITUDE tells "how much" of a certain frequency component is present and the PHASE tells "where" the frequency component is in the image. I understood about the MAGNITUDE part but not able to understand PHASE part. Please explain what it means by "PHASE tells "where" the frequency component is in the image". Thanks


texts: Can a D2D user (a device) transmit simultaneous different signals to a cellular user and a D2D receiver. If we suppose that the D2D device is equipped with SDR technology and can change their frequency channel while transmitting to a cellular user. Moreover will the time slots and power factors with which it is transmitting will be different in this case.


texts:
I'm writing a viterbi decoder and I'm rather confused about the correct process of decoding the final k bits for a kth order convolutional filter.

For the bulk of the message, my process considers the 0-successor and 1-successor of every state, that is, the resulting shift register state found by shifting a 0 in or a 1 in, respectively. Both of these successors share a least-error ancestor found by the dynamic programming step in the Viterbi algorithm.

For the last k bits, we know that the encoder shifted in 0s in order to flush out the payload and reset to the 0 state. I can confirm this is the case because I've also written the encoder. So, when I get to this part of the message, my decoder only considers 0-successors. I think this is all I can do with this information, but I feel like I'm possibly missing some extra info that would help recover the message.

What I'm discovering in reality is that unrecoverable errors are far, far more likely to occur in the final k bits of the message leading up to these 0 flush bits. For example, if my decoded message is 256 bytes long, and I inject errors at a rate of 3%, I find that 15% of the unrecovered errors occur in the last byte.

To some extent, I might expect this to be the case as there aren't additional bits left to help the Viterbi best path converge. I'm using a convergence length of 5k, or 35 decoded bits for a 7th order filter. Should I just expect that without a padding of 35 bits, the end of my message is more likely to receive errors? Or should I keep looking for a bug in my decoder implementation? In practice, is it common to add additional payload bits to help viterbi converge?


texts: phenomenon when window size for STFT is too short

texts:
Hi I was researching on STFT when I encountered this slide online

I don't understand the points, Can anyone explain them? Thanks in advance.


texts: I am trying to generate a sine wave digitally. I am using LMS7200m RF chip from Lime microsystems. I want to generate a sine wave of 100KHz. The DAC on the chip has a sampling rate of 640 MHz with 12bits resolution. How can I do it in C++ or C code reliably? The Rf chips accepts only I and Q values as inputs. I am new to digital signal processing can any one help me understand on how to generate I and Q values for the sine wave in question. Thanks.


texts: How to reconstruct the signal if there were some temporary offsets

texts:
I am not a signal processing expert. I have data like this:

where the red is what I have and the black is what I want to have. Any idea how to reconstruct this? The code for generating the data is here:

close all
Ndata = 1000;
X = nan(Ndata,1);
sigma = 0.5;
Y(1) = 50;
offset = 0;
for i =2:Ndata
    % consumption
    X(i,1) = -0.8*rand;
    if rand<0.01 % supply
        X(i,2) = 50;
    end
    if rand<0.5 % offset change
        if abs(offset)==0
           offset = exprnd(50)*(randi(2)*2-3);
        else
           offset = 0;
        end
    end
    X(i,3) = offset ;
    % noise
    X(i,4) = sigma*randn;
    Y(i) = Y(i-1)+sum(X(i,1:2));
    Z(i) = Y(i)+sum(X(i,3:4));
end
plot(Z,'r')
hold on
plot(Y)

The tricky thing is the fact that there is not only a small white noise, but there are also temporary excesses that are relatively large.

Thank you for any advice.


texts: Algorithm(s) to mix audio signals without clipping

texts:
I'd like to mix two or more PCM audio channels (eg recorded samples) digitally in an acoustically-faithful manner, preferably in near-real-time (meaning little or no peek-ahead).

The physically "correct" way to do this is summing the samples. However when you add two arbitrary samples, the resulting value could be up to twice the maximum value.

For example, if your samples are 16-bit values, the result will be up to 65536*2. This results in clipping.

The naive solution here is to divide by N, where N is the number of channels being mixed. However, this results in each sample being 1/Nth as loud, which is completely unrealistic. In the real world, when two instruments play simultaneously, each instrument does not become half as loud.

From reading around, a common method of mixing is: result = A + B - AB, where A and B are the two normalized samples being mixed, and AB is a term to ensure louder sounds are increasingly "soft-clipped".

However, this introduces a distortion of the signal. Is this level of distortion acceptable in high-quality audio synthesis?

What other methods are there to solve this problem? I'm interested in efficient lesser-quality algorithms as well as less-efficient high-quality algorithms.

I'm asking my question in the context of digital music synthesis, for the purpose of mixing multiple instrument tracks together. The tracks could be synthesised audio, pre-recorded samples, or real-time microphone input.


texts: conversion of 1d signal to 2d signal for CNN operation

texts:
I have a 1D signal data sampled at 25kHz and I want to modify in to 2D signal so as to modify it as the input to ConvolutionalNets as they work only on 2 dimensional data.
What are possible ways to convert the signal from 1D to 2D for such calculations. I have seen some papers but I couldnt find any relevant understanding how the conversion of signal is modified as to apply to CNN operations.
For example if I have 25 data samples in the input and if i rearrange in to a 2d signal as 5X5 sized(inbuilt function in matlab: "reshape").If I am applying 3x3 filter for convolution operation, smallest filter and the output of that 5x5 size input will be 3x3 size.
Can I apply such input to any convnets? Is that a valid type of conversion to apply for CNN ? In such kind of inputs I understand there will be relation in row wise but not column wise.
Any explanation or suggestion would be helpful.

Updated for the comment

I want to convert 1d signal in to 2d image signal. In an image a pixel will have relation with the adjacent 8 pixels. If I am rearranging a 1D signal to 2D signal as quoted in the example , the row 1 will have 5 data samples and next row starts with 6th sample and so on. If I take a sample in row 2, 2nd sample, I have relation first and third sample in the same row but nothing with samples from row 1 and row 3. So If Im converting such way is that appropriate ? In the paper mentioned ( Fig 1) I wanted to know how the 1D signal is converted and applied to CNN link


texts:
I know this is a bit of a fundamental DSP question but I want to finally understand this and I think I need some help to finally put it all together.
Let's say I have an RF complex sampling ADC with a sample rate of 4GSPS.  The device has built in decimation with a minimum decimation of 4X so my maximum capture bandwidth is 1GHz (I know slightly less to allow for aliasing guard bands etc.).
So I would put a 1GHz analog anti alias bandpass filter in front of the ADC to avoid any aliasing, and capture the full 1GHz bandwidth with the ADC. At this point I effectively have 4X oversampling and my data is output to the FPGA at 1000MSPS IQ samples.
The question I have is if I now want to break up this bandwidth in the digital domain into smaller pieces for simultaneous processing of different channels, lets say 100MHz chunks for simplicity.
Is that possible with decimation without breaking Nyquist rules?  My analog anti alias filter is still 1GHz, it would seem like if I decimate by 40x to get the 100MHz bandwidth channels from the original samples that I would now have an effective sample rate of 100MSPS.
Does this work because the original data was taken with a higher sample rate and I am decimating and digital filtering non-aliased data?
I think that is the root of what i was getting at, would those filters need to be analog and limit the BW to half the final sample rate after decimation? The real question ends up being: is it possible to sample a wide-bandwidth (1GHz) and break it up digitally by digital filtering and decimating?


texts:
I'm not sure how to approach this problem.  I will describe what I am trying to do, and what some of my matlab outputs look like.  I'm not trained in signal processing or anything, so please be patient reading.

I have a fourier transform of some signal, and its frequency domain covariance matrix.  For example, say my signal $x[f]$ is $100$ points long.  Its covariance matrix $\Sigma$ is $100 \times 100$.  I noticed that its diagonal is real, and when I plot it, it looks like a PSD.  All the off diagonal elements are complex.

What I wish to do, is take the IFFT of the signal $x[f] \rightarrow x[t]$, and then plot $x[t]$ with error bars from the covariance $\Sigma$.  Is this possible?  After I take the IFFT of $x[f]$, what do I need to do to $\Sigma?$

Thank you for reading!


texts: Oversampling then filtering and downsampling

texts:
Currently I'm working on a project to record ELF (Extremely Low Frequency) signals from a receiver. The highest frequency of interest is about 30 Hz and the system should have a flat frequency response.  I'm looking for a simple solution using an 8 bit ADC and I wonder how the resolution can be improved using oversampling and to alleviate the requirements for the analog low pass filter before the ADC. Without oversampling and a sampling rate of 64 sps a rather steep analog low pass filter is needed to avoid aliasing. With oversampling and a gain of 3 more bits the sampling rate is 2^(2*3) times 64 sps = 4096 sps and the analog filter before the ADC is much easier to realize. The problem is to convert the 4096 sps signal back to 64 sps using a digital low pass filter and downsampling.

From what I read so far, the attenuation needs to be almost 68 dB at 32 Hz for 11 bits but at the same time signals up to about 30 Hz should pass. It seems a very steep digital low pass filter is needed. I tried to do some calculations with filter design software, but the required order is too high to get results. Is this really so hard to realize or am I missing something? Any help would be greatly appreciated.


texts: Choosing the right overlap for a window function

texts:
On Wikipedia there is quite a list of window functions.  However, when you are doing a lot of STFTs, you often have to have some overlap to get better results.

Generally an overlap factor of 2 (i.e., 50% overlap) works fine (I'm doing analysis & resynthesis), but for some windows it isn't good (I currently have this experience with the Blackman window).

I tried my best but I couldn't find a resource that would list the "good" overlap factors for common and less common windows.


texts: How to Show A OFDM Symbol's PSD in Matlab

texts:
I'm trying to suppress the sidelobe of OFDM but having troubles showing its spectrum in Matlab. After IFFT, I get a time-domain signal, but in which way I can get its spectrum(or PSD)? I want to analyze the out band power emission, with rectangle window in the time domain the out of band power curve decrease slowly, so it may interference the adjacent channel user, to avoid that situation we have to apply some algorithms to suppress the out of band power emission caused by the sidelobe of OFDM signal. After apply the algorithm I have to plot some figures to demo the effect of my algorithm, so I need the whole PSD of OFDM signal and I also want the out of band part become more steep. I tried to use FFT and pwelch function of Matlab, but never get what I want. I doubt if I use them in right way. So I'm hoping someone here can help me. I know it's basic, but it really obsess for a long time. I analyze the baseband signal, so it should be after the IFFT operation.The figure that I want is below.

My matlab code:

  N = 64;
  over_sample_factor = 2;
  M = N*over_sample_factor;
  Mod = 16;
  symbol = 1;
  bitlength = N*log2(Mod)*symbol;
  itr_num = 1000;
  fft_len = 2*M;
  signal_freq = zeros(itr_num,fft_len);
  for itr = 1:itr_num
      bit_data = randi([0,1],bitlength,1);
      h = modem.qammod('M', Mod, 'SymbolOrder', 'Binary', 'InputType', 'Bit');
      zp_before = modulate(h,bit_data);
      after_zp = zeros(1,M);
      after_zp(1:N/2) = zp_before(N/2+1:N);
      after_zp(M-N/2+1:M) = zp_before(1:N/2);
      ofdm_symbol = ifft(after_zp);
      signal_freq(itr,:) = abs(fft(ofdm_symbol,fft_len)).^2;
  end
  PSD_mean = mean(signal_freq,1);
  plot(fftshift(10*log10(PSD_mean)));

The figure is like this


texts: How to convolve an arbitrary signal with a causal decaying exponential?

texts:
I need to simplify the following convolution

$$x(t)\star [e^{-2t} u(t)]$$

where $u(t)$ is the unit step function. I'm very confused with this. Using the definition of convolution of continuous-time signals, I obtained

$$\int\limits_{-\infty}^t {x(\tau)e^{-2(t-\tau)}u(t-\tau)}{d\tau}
$$

However, I have no idea what to do at this step. I'm stuck. Is there any property / theorem / trick that can help me simplify this?


texts:
I'm trying to find an intuition behind auto-encoder using an analogy with filter banks.
I can comprehend the encoder and analysis filters in a filter bank as extracting features from the input signal i.e. they represent the signal in more and more complicated features in deeper layers, and the original image is like a combination of those feature "basis". However I'm trying to get an intuition behind the reconstruction of the signal. The decoder and the synthesis filters only have a lower dimensional representation of these signals, usually which is very sparse, then how do the reconstruct intuitively?
Are the analysis and synthesis in a filter bank inverses of each other? In order to achieve the best reconstruction is there any relationship when choosing the synthesis filters? I don't think they can be inverses because they operate on different data, i.e. the synthesis filters get 0 padding which wasn't there originally, so perhaps an inverse would not provide the best reconstruction or will it?
Similarly, are the transposed convolutions in the decoder part of the auto-encoder inverse of the convolutions in the encoder? Some papers use the decoder as the literal transpose of the filters in the encoder. I don't think they both can be inverse operations right just like transpose of a matrix is not inverse unless in some particular scenarios? If the weights are not shared between the encoder and decoder, what exactly is the decoder trying to do? Is it really trying to learn the inverse of the encoder? Because if that's the case then there shouldn't be a need to have the decoder learn anything, we could simply always use decoder as the inverse of the encoder throughout the training and learning?
Is there a similarity in this sense between filter bank and auto-encoder?


texts: What is meant by LFM (Linear Frequency Modulation) waveform?

texts:
I m a newbie in radar experiment on which currently i m working, and I m trying to understand the transmitted waveform which is LFM (Linear Frequency Modulation), but unable to interlink between the following parameters of transmitted signal waveform ,can someone explain to me how do these are connected to LFM?

Center frequency  : 150kHz
Sampling Frequency: 4000kHz
Bandwidth         : 50kHz
Pulse length      : 500microseonds


texts: How does sampling and subsequent resampling affect aliased harmonics?

texts:
If I sample a signal with a harmonic higher than half my sample frequency I can predict the aliased frequency easily as

$$|F_h-kF_s|$$

But what if after sampling the signal, I re-sample at a lower rate (decimate)?

Does the lower sampling frequency simply trump the original higher frequency or will there be mixing or modulation effects due to the two sampling frequencies that change or create new aliased harmonics?

If so, is there a way to predict what the aliased harmonics will be?

As a simple example assume I sample a signal containing a $7$ kHz harmonic at $3$ kHz, then decimate by taking every third sample, so effectively re-sampling at $1$ kHz.


texts:
What is the RFFT and iRFFT complexity for 1D and 2D?
Here is some starting point.
1D RFFT

The FFT is $\mathcal{O}(N\log N)$
RFFT computes $N//2 + 1$ frequencies instead of $N$.

Thus, RFFT is $\mathcal{O}([N//2 + 1]\log N)$ = $\mathcal{O}\left (\frac{N}{2}\log N \right )$

Question: What about iRFFT? Is it $\mathcal{O}(N\log N)$ or $\mathcal{O}\left (\frac{N}{2}\log N \right )$ ?

2D RFFT

The 2D-FFT is $\mathcal{O}(2N^2\log N)$
RFFT computes $(N//2 + 1)N$ frequencies instead of $N^2$.

Thus, RFFT is $\mathcal{O}([N//2 + 1]N\log N + N^2 \log N)$ = $\mathcal{O}\left (\left [ \frac{3}{2} N + 1\right ] N \log N \right )$ = $\mathcal{O}\left (\frac{3}{2} N^2 \log N \right )$

Question: What about iRFFT? Is it $\mathcal{O}(2N^2\log N)$ or $\mathcal{O}\left (\frac{3}{2} N^2 \log N \right )$?

Am I getting anything wrong?


texts:
I'm reading the book "Digital Image Processing" by Gonzalez and Woods and I'm wondering how their definition of zero-phase-shift filter is equivalent to the one given here.

"Digital Image Processing" gives the following definition:

The DFT of $f(u, v)$ in Polar coordinates is:

$$F(u, v) = \sum_{x=0}^{M-1}\sum_{y=0}^{N-1} e^{-j 2\pi (\frac{ux}{M} + \frac{vy}{N})} = |F(u, v)|e^{j\phi(u, v)}$$

where $\phi(u, v) = \arctan\left(\frac{I(u, v)}{R(u, v)}\right)$ and $I, R$ is the imaginary and real part of $F(u,v)$. Then filtering means multiplying the filter $H(u, v)$ with $F(u, v)$ and then taking the inverse DFT i.e. $F^{-1}\left(H(u, v)R(u, v) + j H(u, v)I(u, v)\right)$.

"The phase angle is not altered by filtering in the manner just described because $H(u, v)$ cancels out when the ratio of the imaginary and real parts is formed. Filters that affect the real and imaginary parts equally, and thus have no effect on the phase, are appropriately called zero-phase-shift filters.".

If I understand it correctly, they are saying $\phi(u, v) = \arctan\left(\frac{\text{Im}(H(u,v))I(u, v)}{\text{Re}(H(u,v))R(u, v)}\right) = \arctan\left(\frac{I(u, v)}{R(u, v)}\right)$.

But the other definition looks different. They focus on the imaginary part of $H(u, v)$ being zero. This makes more sense to me, because then $\phi(u, v) = 0$. Right? But how does this affect "the real and imaginary parts equally"?


texts: Looking for good bood suggestion relate dsp and antenna thoery

texts:
I am a new student in dsp field and looking for a good book to practice at home.

Please provide good books.

Thank you.


texts: Relation between sawtooth Fourier coefficients and its DFT

texts:
I'm having some trouble with understanding the DFT of a sawtooth single period signal and its relation with sawtooth Fourier coefficients.

Let's say I have a signal
$$
s(t) = \frac{At}{T} - \frac{A}{2} \qquad t\in[0,T)
$$

If I plot its DFT real and imaginary components I have something like this

Where $A=10$, $T=1$, $N=50$ (number of samples) and the DFT is normalized by $\tau = T/N$ and $1/N/\tau$ for the inverse.

Now, the imaginary part as far as I can understand comes from the Fourier coefficients
$$
c_n = \frac{iA}{2\pi N}
$$
These coefficients (see e.g. Brigham, section 5.2) equal the fourier transform $S(t)$ scaled by $\frac{1}{T}$ and sampled at $n/T$. Given we have a discrete signal, we need to account for aliasing in the frequency space, which can be represented by repeating $c_n$ at $N$ intervals and summing with themselves, something like the following? probably missing a constant
$$
Im(S_n) = \sum_{m=-\infty}^{\infty} c_n\delta(n-mN) \qquad n \in [0,N-1]
$$

What about the non zero real part? As an afterthought I can rationalize it because the sawtooth is non zero at $t=0$, so it must have a non zero mean DFT. But if I take a sawtooth centered in $-T/2,T/2$ it's still there even if the signal crosses the origin.

I'd like to find a better relation that gets me from the Fourier series to the DFT, accounting both for the imaginary and the real part. Any pointer to good resources to help understand the whole thing would be more than welcome.

Just looking at the DFT, the real part, with my normalization should be:
$$
-\frac{A}{2}\tau\sum_{n=-\infty}^{\infty}\delta\left(f-\frac{n}{T}\right)
$$
How do I justify this result? especially the multiplication factor?

(sorry if the math is a bit flaky, I'm still trying to understand everything about DFT, Dirac combs and most of all normalization constants, every single source has them different...)


texts: Explanation of FIR and IIR Filters with Simple Words

texts:
It seems I am unable to understand how these filters work and what is the point of using them, no matter how many tutorials I watched.

Could someone please explain, in simple words, what FIR and IIR filters do?


texts: What is the easiest learning curve to develop a Mark Reader software?

texts:
I am a C# programmer and have a project to create a Mark Reader software. Bubble sheets that have been filled in by students are scanned using a general scanner (as opposed to a special scanner for bubble sheet). The resulting images must be processed using the Mark Reader software in question.

I have a bit knowledge in DSP. Shortly speaking, what is the easiest learning curve to develop the software? Book recommendation is also welcome.


texts: Estimating a Signal Given a Noisy Measurement of the Signal and Its Derivative (Denoising)

texts:
I have a signal and its derivative simultaneously measured, both including additive noise. The measurement is completed before the analysis, so it can be looked ahead. Now I want to reconstruct a less noisy version of the signal. I'm looking for pointers to algorithms I should look into.

Kalman filter seems to be on the right track, but the implementations I see so far are trying to estimate based on previous measurements only, while I probably should use both previous and coming measurements for optimal results at each point.

Ideas?


texts: Smoothing signal / detecting bumps in a data stream

texts:
(EDIT: This question follows from Extracting Binary Magnetic-Strip Card Data from raw WAV)

Here is my signal (top line) and a basic IIR filter applied (bottom line)

(EDIT: my task is to break the signal into binary 0 (frequency F) and binary 1 (frequency 2F) -- that's why it is called F2F.  So I need to process it in such a way that guarantees no false peaks.  While the screenshot makes it look trivial, there is a potential problem of getting a double peak, and also of getting false positives in the trough between real peaks.)

My question is, what methods are available for smoothing this signal?  Is IIR my best bet?

I can see at least three possibilities:

IIR y[n] = 0.9*y[n-1] + 0.1*x[n]  where y[x] =0 when x < 0
Moving / windowed average -- place a Bell curve with area 1.0 over the surrounding say w=10 samples each side and integrate bellSmooth(x) = integral[x-w,x+w] { bell(k).samp(k) }dk
Determine the expected frequency and FFT / remove higher order bins / reverse FFT

I may have answered my own question, but probably this is incomplete and I'm sure I am using the wrong terminology.  Also I can't really predict the pros and cons.   The last method is less attractive as it requires knowledge of the basic signal frequency.  But then so does the second method; I need to choose an appropriate window length.

Are there any other methods?


texts: Low pass filter output showing instability

texts:
I am fairly new to DSP, and I'm trying to implement the Pan Tompkins algorithm for QRS detection of ECG signals in MATLAB. The first stage of the algorithm consists of a second-order low-pass filter. The difference equation of the filter is

$$
y(n)=2y(n-1)-y(n-2)+x(n)-2x(n-6)+x(n-12)
$$

Here is the MATLAB code I have written.

load 108m.mat; sig_108 = val(1,:); clear val;
lowPassResult = lowPassFilter(sig_108);

figure;
subplot(2,1,1)
plot(lowPassResult);
title('Lowpass')
subplot(2,1,2)
plot(sig_108);
title('Original')

function y = lowPassFilter(x)
    %y(n)=2y(n-1)-y(n-2)+x(n)-2x(n-6)+x(n-12)
    y = zeros(1,21600);
    y(12) = 0;
    y(11) = 0;
    for n = 13:length(x)
        y(n) = 2*y(n-1)-y(n-2)+x(n)-2*x(n-6)+x(n-12);
    end
end

I am using ECG record 108 from the MIT-BIH arrhythmia database available at PhysioBank ATM. One problem I am having is that the output of the lowpass filter appears to grow without bound towards $+\infty$ as can be seen in the figure above. Where have I gone wrong in my code that could be causing this issue?


texts: My teacher said to me that Instantaneous Angle means Argument of Complex Function thats all. Someone please tell me what significance does this "instantaneous" Word has here ?


texts: Is there a way to construct FIR analog low-pass filter? (By this, I believe I am asking for no-pole analog low-pass filter, in terms of Laplace transform.)


texts: Nyquist plot interpretation when curve hits the origin

texts:
I'm a bit confused about the interpretation of the Nyquist plot when the origin is part of the plot. In this case I'm not even considering closed loops, I'm just looking at the Nyquist Plot of a given transfer function and trying to relate it to the corresponding poles and zeros.

For example, if I take the Nyquist Plot of $\frac{1}{s}$:

This system is marginally stable, and the origin is neither inside nor outside the loop, it is in the curve.

Now if I take Nyquist Plot of $\frac{1}{s+1}$:

Again, the origin is not outside the curve, but now the system is stable.

And, lastly, if I take an unstable system, like $\frac{1}{s-1}$:

The origin is, again, part of the curve.

In my understanding, the origin of the Nyquist Plot of $\frac{1}{s+1}$ should be outside the curve, since it has no poles and zeros in the right half plane. Also, the origin of the Nyquist Plot of $\frac{1}{s-1}$ should be surrounded by a counter-clockwise loop, since it has 1 pole and 0 zeros in the right half plane.

Shouldn't the Nyquist Plot of a given transfer function tell the difference between poles and zeros just by looking at how many times (and which direction) the origin is surrounded by plot loops?


texts: If a function has slowly varying amplitude and phase what is the bandwidth of the amplitude and phase individually?

texts:
I am dealing with time histories that are narrow band with slowly varying amplitude and phase. I wish to extract the amplitude and phase but there is noise on the time history. An example of the type of time history I am dealing with is as follows:

This is a decaying oscillation but does not have an exponential decay. A short part of the time history shows that there is noise on the decaying amplitude:

The modulus of the Fourier transform of the time history shows a clear peak:

I have drawn lines at $f_1$ and $f_2$ that contain the bandwidth of the time history I am interested in. However, there is noise at a frequency of about 230 Hz that spoils the data.
If I extract the amplitude and frequency the noise can be clearly seen.

I wish to remove the nose from the amplitude and frequency time histories I have estimated. Wikipedia gives many examples of smoothing methods such as low pass filter, Savitzky-Golay and moving average. Each of these methods requires a knowledge of the bandwidth of the signal being smoothed. My question is: what is this bandwidth?
My intuition is that the bandwidth is equal to 0 to $f_2-f_1%$ Hz or a smoothing time interval of $\frac{1}{f_2-f_1}%$. Is this correct? Is there a rigorous proof of this or are there references you can point me towards?
I think the problem is different to finding the amplitude of an amplitude modulated frequency because my frequency is also time varying.
Thanks for any help.


texts: Blind source separation on images - with a known source

texts:
A camera outputs an image $A$ that is a linear combination of two uncorrelated images $B$ and $C$:

$$A = B + kC$$

Images $A$ and $C$ are known; the unknowns here are $B$ and a scalar coefficient $k$ that varies the intensity of $C$.

Is there a way to retrieve $B$? I could use $A$ and $C$ as separate sources and apply the ICA algorithm, but I think that it's an overkill in this case, and there should be something simpler. Any suggestion? (I'm using OpenCV to process the images, but I can implement non-OpenCV algorithms if I need to.)


texts: Deciding if the received packet is valid or not valid

texts:
I have sampled RF wireless packets as IQ samples file, those RF packets has patterns like preamable|syncword|data|crc , I'm doing cross correlation with my constant syncword (after I did the demodulator in order to convert the IQ samples ) and I see in matlab the peak of correlation where there's sync word in my packets, Im now trying to decode those RF packets which I want to decide if the packets that I received are valid or not valid.
what should I do in order to decide if the packet is valid or not?? maybe it's related to CRC but how can I decide if the packet is valid or not? what I have done and arrived to this stage : I did cross correlation with my syncword and I see the peak of corelation of syncwords !
Anyone could help me by the flowchart what should I do in order to decide if my received packet are valid or not valid? thanks
CRC/SYNCWORDS in the transmitter are constants ..


texts: Band edge filter for raised cosine impulse

texts:
I'm hoping to use an implementation of a frequency-locked loop for rough frequency synchronization in a PSK31 demodulator. The approach is to define a filter that is the derivative of the matched filter in the frequency domain. A couple slides from fred harris illustrate the idea:

I've found implementations for the typical case of a root-raised-cosine pulse shaping filter. However, PSK31 is a bit "special" in that it uses a raised cosine impulse. AKA, the Hann function:

$$ h(t) = {1 \over 2}\, (1+\cos(\pi t))\, \Pi(t/2) $$

Where $\Pi$ is the rectangle function.

How can I calculate the appropriate band-edge filter in this case?


texts: Intuitive meaning of complex time domain signal representation (quadrature sampling)

texts:
I saw What is the meaning of imaginary values in the time domain?.

My question is about I Q samples as recorded from an SDR receiver (e.g. RTL-SDR).

My understanding is that we can call these I Q samples a complex time domain signal and these can have both non-zero real and imaginary components (the question I linked above is talking about zero imaginary part in the time domain).

Is there an intuitive explanation of what is the meaning of this type of complex time domain signal?

I can think of taking a simple case of an input sine wave of a frequency below the sampling frequency which would make the IQ signal's phasor rotate with the frequency equal to the difference between the input sine wave's frequency and the sampling frequency thus allowing to unambiguously reconstruct frequencies up to (but not including) the sampling frequency from the IQ samples.

But I'm struggling to extend this special case to an input composed of several sinusoids of different frequencies.


texts: Spectral density of modulated noise process

texts:
Suppose we have a noise process $F(t)$.
Suppose $F$ is stationary and we know its autocorrelation function
$$\rho_F(\tau) \equiv \langle F(0)F(\tau) \rangle$$
where the average is over realizations of the process.
We then construct a new noise process $G(t) \equiv \cos(\Omega t)F(t)$.
Obviously, $G$ is not stationary.
The autocorrelation function of $G$ is
\begin{align}
\langle G(t_1) G(t_2)\rangle
&= \langle F(t_1) \cos(\Omega t_1)F(t_2) \cos(\Omega t_2) \rangle \\
&= \langle F(t_1) F(t_2) \rangle \cos(\Omega t_1) \cos(\Omega t_2) \\
&= \rho_F(t_2 - t_1) \cos(\Omega t_1) \cos(\Omega t_2) \, .\\
\end{align}

The spectral density of $F$ can be usefully defined as
$$S_F(\omega) = \int \frac{d\omega}{2\pi} \rho_F(\tau) e^{-i \omega \tau} \, . $$
However, since $\langle G(t_1) G(t_2) \rangle$ depends on $t_1$ and $t_2$ and not just their difference, we cannot easily define $S_G(\omega)$.
This is distressing because for deterministic signals we expect modulation by a sinusoid to just shift the spectral density, so we expect there should be a simple formula for $S_G(\omega)$.

What, if any, is a meaningful way to construct the spectral density of $G$?


texts: Threshold Estimation from a Linear FM Pulse Position Modulated signal

texts:
I have a signal chain to detect multi channel LFM (Linear Fm - Hann window being applied as shaping filter for the pulses) pulses as shown in picture below. (Numbers in the picture is in Hz., for instance start frequency of channel-1 LFM is 39400 Hz).

As you may see I have one linear chirp signals per channel sweeping 300 Hz in 60 msec with different start frequencies (meaning, if you look at the block diagram channel-1's LFM pulse starts at 39400 and sweeps all the way up to 39700 then channel 2's LFM pulse LFM pulse starts  at 39800 and sweeps all the way up to 40100. This is a multichannel system so 16 channels can be active anytime simultaneously.
This signal is going through a Analog to Digital Converter Circuit with a sampling rate of 500 kHz. I have no control on this circuit. I can only play with back-end DSP algorithms which are depicted with blue block in the picture.
So ADC board mentioned above, the does a complex filtering and down samples the signal by 4 then another complex filter and down samples the signal by 2. So when I receive the signal from ADC board it comes through UDP with 8192 samples in each transfer with the sampling rate of 62.5 KHz (8192 real and 8192 imaginary symbols).
As you can see I am applying band-pass filter for each channel to achieve better channel separation. and this block is followed by another down-sampler. After the match filters (which are also designed with considering this decimation by 16), I am getting peaks at where match filter "matches". Now I need to find "the potential peak locations and their amplitudes" to feed it to my peak detector algorithm. Eventually I should be able to calculate the time difference between 2 LFM pulses on each single channel to decode the information carried by these pulses.(Pulse Position Modulation).
I want to be able to estimate the "noise level" so that I can place a threshold to determine potential peak candidates. I should have used another word for the noise floor since this is using purely time domain signal. Anyways, I want to calculate the "floor level" and then estimate a threshold that will be half way up through from the estimated "floor" (With the measurements we are confident that half way up would be Enough).
So far I am working on a method that does the following:
th = mean(movmedian(20*log10(inputsignal),3));
I was hoping to get rid of peaks with moving average filter and mean value of this remaining signal should give me a coarse value close to "floor level". I am not sure about how big the median window should be?
Can someone comment on a new method or the one I am applying now?
Thank you in advance!


texts: Working backwards from PSD to possible signal

texts:
I have been trying to reconstruct a random signal from its PSD and am running into trouble. I know that many different signals in the time or spatial domains can result in the same PSD-- I am interested in any of these signals. I am not the greatest at signal processing, so I've gotten a bit confused. I am reading this paper. They state:

  Therefore, PSD is approximated as
  $$
S(k)=\frac{2\sigma^2L_c}{(1+k^2L_c^2)^{0.5+\alpha}}
$$
  where
  $$
k=2\pi f, f=i/(N\Delta Z), 0\lt i \lt N/2
$$
  Lastly, N is the number of points along the line.

  ...

  With the magnitude information provided by S(k), we can reconstruct random line edges by applying a random phase to each frequency component of the PSD to form a unique signal in the frequency domain. A line edge with roughness can be simulated by doing an inverse Fourier transform of this signal. Random lines are distinguished through applied random phases.

I have computed values for S(k), per the equation given. I then tried taking the square root of those values, multiplied them by the number of values, and subtracted their mean. Finally, I generated random phase information before doing an inverse FFT. Even after subtracting the mean I still see a very large spike at the beginning and end of the signal.

Does anybody have any suggestions as to what I'm doing wrong?

For those that speak Python, here is the code I've written so far:

import numpy as np
from matplotlib import pyplot as plt

if __name__ == '__main__':
    SIGMA = 4.0
    ALPHA = 0.2
    Lc = 25

    L = 2000.0
    N = 1000

    dz = L/N

    i = np.linspace(0,N/2,N)
    f = i/(N*dz)
    k = 2*np.pi*f
    PSD = lambda x: (2*3*SIGMA**2*Lc)/(1+(x*Lc)**2)**(0.5+ALPHA)

    window = 25.0/46-21.0/46*np.cos((2*np.pi)/(N-1)*np.array(range(N)))
    magnitude = N*np.sqrt(PSD(k))
    phase = 1j*np.random.randint(0,6,N)+np.random.randn(N)
    FFT = magnitude-np.mean(magnitude)+phase
    FFT = window * FFT
    FFT = np.concatenate((FFT[::-1][0:-1],FFT))

    fig = plt.figure()
    ax = fig.add_subplot(111)

    ax.plot(np.fft.ifft(FFT))

    plt.show()

Edit: I forgot about Gibbs phenomenon! As per Bruce Zenone's suggestion, I applied a Hamming window and mirrored my data. I added plots and updated my Python code. I am still seeing some effects of Gibbs phenomenon. Truncating the affected parts, the signal does not seem to be varying with sigma, alpha, and Lc.

Here is my resulting signal:

And a figure from the paper showing the effect of the different parameters:


texts: Detecting which digital number an image represents

texts:
I photograph the dashboard of my car every day to register the mileage for work.
I use OpenCv to detect the region of the digital displayed mileage, crop to the numbers and match them. The method I use now to determine what numbers are displayed is working for (lets say) 95% of the cases. Finding the numbers is working fine (resize, rotate, correct perspective). Matching the numbers I find still a bit tricky.

The method I use now is determining the center location of each part/stripe that makes up the digit and match the difference of the center to predefined values. To find the center I use blob detection of open cv.

is then converted to a binary image

and blob detection is used to find the center of each blob (marked with red)

The center location of each blob is compared (squared distance) to the center location of a predefined set of locations. And the least square distance is related to what stripe of the digit it is. Then the combination of the stripes found is compared to a dictionary containing the stripe locations for each digit.

The problem I have now is that, although the images are normalized, sometimes more blobs are detected because the the threshold for making the binary image leave out pixels unexpectedly in the center of a blob. And blob detection sees that are separate blobs.

Is there a better approach to determine what digit it is using the particular display of the digits shown above? 'Better' would be: more robust, simpelere, less logic/code required, more universal applicable. I don't want to pull in a whole OCR library for this, if possible.


texts: What is the effect of sampling time, sampling frequency and SNR on latency?

texts:
I believe that if we have two or more external sound cards, hence clock would vary for each device.

How will this affect the latency?
Will it increase or decrease?
How will the factors like

sampling rate,
sampling time,
SNR and
voltage reference impact latency?


texts:
Studying for an exam I'll have I ran into a question about computing the localization of a Gaussian Filter in 1D, but sadly I have no idea of what such localization is. Can you give me a hint?

Thanks


texts:
I'm new to DSP and I would like to make an SDR receiver. For this I want to shift bins in the frequency domain in order to bring to 0Hz any portion of the base band spectrum (I know that there are also other methods). I got this idea from SDR for the masses page 4 figure 7

I use a stereo 192KHz 24bit audio codec connected to a dual channel audio generator with 90 degrees phase shift between channels. I use the complex FFT/IFFT from Paul Bourke FFT with 8192 samples length and no window.

If I do FFT and then IFFT I can't notice any signal degradation (I tested it with a single tone and with music). However, if I try to rotate to left with 100 the bins in the first half of I and Q array and to right the second half of array then I can hear a cyclic noise. The single tone output in this case has some distortions at both ends (second picture), as opposed to the first picture where I did not shift the bins.

Also if I shift the bins in the other direction the single tone frequency also seems increase, but not so much(see bellow)! What am I doing wrong? Thanks

void rotate_left(SAMPLE_DATA_TYPE *source, SAMPLE_DATA_TYPE *dest, unsigned short
    int x, unsigned short int len) {
        memcpy(dest, source + x, (len - x) * sizeof(SAMPLE_DATA_TYPE));
        memcpy(dest + len - x, source, x * sizeof(SAMPLE_DATA_TYPE));
    }

void rotate_right(SAMPLE_DATA_TYPE *source, SAMPLE_DATA_TYPE *dest, unsigned short int x, unsigned short int len) {
            memcpy(dest, source + len - x, x * sizeof(SAMPLE_DATA_TYPE));
            memcpy(dest + x, source, (len - x) * sizeof(SAMPLE_DATA_TYPE));
        }

main() {

................

#define SAMPLES                 8192
#define FFT_SIZE                13
#define SAMPLE_DATA_TYPE        double
shift = 100;

while (1) {

    //Get I and Q from interleaved audio stream
    ..........

    FFT(1, FFT_SIZE, I, Q);
    //I and Q have the DC bin at first and last position. At the center is the highest frequency bin

    rotate_right(I, temp_buffer_I, shift, SAMPLES / 2);
    rotate_left(I + SAMPLES / 2, temp_buffer_I + SAMPLES / 2, shift, SAMPLES / 2);
    rotate_right(Q, temp_buffer_Q, shift, SAMPLES / 2);
    rotate_left(Q + SAMPLES / 2, temp_buffer_Q + SAMPLES / 2, shift, SAMPLES / 2);

    memcpy(I, temp_buffer_I, SAMPLES * sizeof(SAMPLE_DATA_TYPE));
    memcpy(Q, temp_buffer_Q, SAMPLES * sizeof(SAMPLE_DATA_TYPE));

    FFT(0, FFT_SIZE, I, Q);

    //Output I and Q
    ..........

}

short FFT(short int dir,long m,SAMPLE_DATA_TYPE  *x,SAMPLE_DATA_TYPE *y)
{
   long  n,i,i1,j,k,i2,l,l1,l2;
   SAMPLE_DATA_TYPE c1,c2,tx,ty,t1,t2,u1,u2,z;

   /* Calculate the number of points */
   n = 1;
   for (i=0;i<m;i++)
      n *= 2;

   /* Do the bit reversal */
   i2 = n >> 1;
   j = 0;
   for (i=0;i<n-1;i++) {
      if (i < j) {
         tx = x[i];
         ty = y[i];
         x[i] = x[j];
         y[i] = y[j];
         x[j] = tx;
         y[j] = ty;
      }
      k = i2;
      while (k <= j) {
         j -= k;
         k >>= 1;
      }
      j += k;
   }

   /* Compute the FFT */
   c1 = -1.0;
   c2 = 0.0;
   l2 = 1;
   for (l=0;l<m;l++) {
      l1 = l2;
      l2 <<= 1;
      u1 = 1.0;
      u2 = 0.0;
      for (j=0;j<l1;j++) {
         for (i=j;i<n;i+=l2) {
            i1 = i + l1;
            t1 = u1 * x[i1] - u2 * y[i1];
            t2 = u1 * y[i1] + u2 * x[i1];
            x[i1] = x[i] - t1;
            y[i1] = y[i] - t2;
            x[i] += t1;
            y[i] += t2;
         }
         z =  u1 * c1 - u2 * c2;
         u2 = u1 * c2 + u2 * c1;
         u1 = z;
      }
      c2 = sqrt((1.0 - c1) / 2.0);
      if (dir == 1)
         c2 = -c2;
      c1 = sqrt((1.0 + c1) / 2.0);
   }

   /* Scaling for forward transform */
   if (dir == 1) {
      for (i=0;i<n;i++) {
         x[i] /= n;
         y[i] /= n;
      }
   }

   return(1);
}


texts: Why Does FISTA Algorithm Not Work for Signed Signals?

texts:
Using the FISTA Algorithm for compressive sensing from https://github.com/tiepvupsu/FISTA, I created the matlab example below.

I created 2 sparse signals x_signed and x_pos, where the latter only contains positive values.
Using opts.pos = false as input to the FISTA, I assume, that the algorithm also works for signed input signals.
As visible, the FISTA algorithm is able to perfectly reconstruct x_pos, but not x_signed.

s = rng(5);
close all;
nrbins = 5;
N = 120;
m = 24;
%create input signals
x_pos = [round(10*rand(nrbins,1));zeros(N-nrbins,1)];
x_signed = [10*randn(nrbins,1);zeros(N-nrbins,1)];

%random shuffle
x_signed = x_signed(randperm(size(x_signed,1)));
x_pos = x_pos(randperm(size(x_pos,1)));

%Compressive Measurements on positive signal
A = 0.5*(sign(randn(m,N))+ones(m,N));
y = A*x_pos;

opts.pos = true;
opts.lambda = 0.01;
opts.tol = 1e-14;
opts.max_iter = 1000;
xrec_pos = fista_lasso(y,A,[],opts);

%Compressive Measurements on signed signal

A = 0.5*(sign(randn(m,N))+ones(m,N));
ysigned = A*x_signed;

opts.pos = false;
opts.lambda = 0.01;
opts.tol = 1e-14;
opts.max_iter = 1000;
xrec_signed = fista_lasso(ysigned,A,[],opts);

%plot positive
figure;
scatter(1:size(x_pos,1),x_pos);
hold on;
plot(xrec_pos);

legend('x','xrec\_pos');

%plot signed
figure;
scatter(1:size(x_signed,1),x_signed);
hold on;
plot(xrec_signed);

legend('x','xrec\_signed');

How can I reconstruct a signed signal using the FISTA algorithm?

For what reason does the option "pos" exist? Does it have performance implications?


texts: Is $\arctan(t)$ an energy or power signal?

texts:
Is the function $\arctan(t)$  an energy or power signal?
I've tried working out the integrals for energy and power but I get stuck on both.


texts: How to Linearly Combine Two Unbiased Estimators of One Parameter without Knowledge of Their CoVariance?

texts:
I have two unbiased estimators of one parameter, $\tau$. The first estimator, $r_1$, is the better estimator with lower variance than the second estimator, $r_2$.

I also have:
$ \mathbb{E} \left[ {r}_{1} \right] = 0.8 \tau $ , and $ \mathbb{E} \left[ {r}_{2} \right] = 0.2 \tau $.

I want to use both estimators to achieve lower variance. Is there a systematic way to find a linear combination of $r_1$ and $r_2$ that results in a lower variance and unbiased estimate, without a knowledge of the variance and covariance of the estimators?


texts:
Assume FFT of $x[n]$ is $X(k)$ and I would like to find $X(-k)$. I have implemented two techniques:

I am calculating the X(k) = fftshift(fft(x)) and then find X(-k) to be fliplr(X(k)) ;
I write the for loop to correspond each symbol of in $X(k)$ to $X(N-K+1)$ where $N$ is the number of points used for FFT.

I get kind of similar answer but in the second approach I am mapping $X(N)$ to $X(1)$ which is DC, which should not be the case. Therefore I am suspecting I am doing something wrong.

Can anyone tell me what is the correct way to implement $X(-k)$ when you have the $X(k)$ ?


texts: Computing error between filtered and noisy time series signal

texts:
I have a live signal that needs filtering, and I'm trying to compare a Kalman filter with a simple exponential filter. For this I generated an artificial ground truth signal, add some noise and filter it with both methods.

The problem is that, while the result from the Kalman filter 'looks' better, I can't quantify how much better it is because its output lags behind the signal. When I simply compute the error between a noisy sample and a filtered sample at time T, the difference is always higher than between the noisy sample and the ground truth sample at the same time.

Is there a way of correctly computing the error in these cases? I thought about computing a least squares fit between the filtered and ground truth signals, but that feels like 'cheating' a bit, as it looks at the best possible case and also doesn't account for the filter lag.

Edit:
Some graphs to give a better picture of what I'm doing

I've plotted the outputs up to t=140. Both the exponential filter (red) and kalman filter (green) are lagging behind the current noisy sample (blue). If I try to compute the squared error by comparing the ground truth value at t=140 with kalman or exponential values at t=140, the error will be large, and at times even larger than the error between the ground truth and the noisy signal (second picture top).

This makes the exponential filter look better, just because it has less lag than the kalman filter.
If I try to account for lag, i.e. compute the error between the kalman filter at t=140 and ground truth at the closest t ~ 137-ish, then the squared error graphs look better (second picture bottom), but I'm not sure if this is the correct thing to do.

Edit 2
I now realize that these are not really time series, but rather processed as a 2D signal. So 't=' should really be 'x=', with the value being y. This explains why the sequences are not synchronized in time. But I think the question is still valid.


texts: Why would FFT interpolation using zero-padding undershoot the true frequency of a single tone sinusoid?

texts:
Why would FFT interpolation by zero-padding or using the Chirp Z-Transform produce a maximum at a bin that corresponds to a frequency less than the input frequency of a single tone sinusoid?

I am attempting to precisely determine the frequency of a single tone sinusoid from a data set which captures just 1-2 wavelengths of the tone. To do this I have tried to use  zero-padding interpolation and the Chirp Z-Transform. When I perform my interpolation the maximum in the FFT or CZT falls into a bin that does not correspond to the frequency of the input sinusoid and instead always undershoots. Further, a convolution between the input sinusoid and the complex exponential of the same frequency is smaller in magnitude than the convolution between the input and the complex exponential of the frequency corresponding to the maximum FFT bin. This obviously goes against the properties of sinusoids.

I have tried to diagnose this issue by windowing, using different interpolation factors, using different length and frequency input signals, and using other fft codes, to no success.

Below is the matlab code that performs the interpolation of the fft of an input sinusoid, plots the fft and czt, plots the input sinusoid along with the sinusoid corresponding to the max bin, and finally performs the convolution between the input and the sinusoids of the FFT max bin and the input frequency.

%create input sinusoid
M = 100;
m = linspace(-pi,pi,M);
f = 2;
x = sin(f*m);

%set interpolation factor and run zeropadded fft
N = M*100;
FFT = fft(x,N);

%define czt parameters
r1 = 0;
r2 = .1;
a = exp(2j*pi*(r1));
w = exp(-2j*pi*(r2-r1)/N);

CZT = czt(x,N,w,a);

%Plot interpolated FFT and CZT
figure, plot(abs(FFT))
hold on
plot(abs((CZT)))
legend('fft','czt')

% Find index of maximum bin for FFT and CZT
[~,indexFFT] = max(abs(FFT(1:N/2)));
[~,indexCZT] = max(abs(CZT(1:N/2)));

% Create sinusoids from the extracted bins to compare against input
% sinusoid
a = exp(-2j*pi*(indexFFT-1)*(linspace(1,M,N))/(N));
b = real(a);

a2 = exp(-2j*pi*(indexCZT-1)*(linspace(1,M,N))/(N/(r2-r1)));
b2 = real(a2);

a3 = exp(-2j*pi*(f)*(linspace(1,M,M))/(M));
b3 = real(a3);

figure,plot(linspace(1,M,N),b,'-.')
hold on
plot(linspace(1,M,N),b2,'--')
plot(linspace(1,M,M),b3,':k')
hold off

legend('FFT','CZT','input sinusoid')
title('the sinusoid associated with the maximum FFT bin')

% Perform convolution between zeropadded input and the sinusoid from the
% fft max bin, and compare against the convolution with the true sinusoid with true frequency

x_Padded = zeros(1,N);
x_Padded(1:M) = x;

X = 0;X2 = 0;

for n=1:N
    X =  X + x_Padded(n)*exp(-2j*pi*(f)*(n-1)/(M));
    X2 =  X2 + x_Padded(n)*exp(-2j*pi*(indexFFT-1)*(n-1)/(N));
end

X2_mag = abs(X2);
X_mag = abs(X);


texts: How do we calculate bandwidth from signal formula?

texts:
How can we calculate the bandwidth of distinct signals like this?
I guess the right answer is 90 - 10 = 80 Hz. But I have a hard time understanding why is this calculated in that way. I thought the bandwidth is calculated by the diversity of signals and since there aren't in-between frequencies like 11, 12... in the image there are 5 signals, shouldn't the bandwidth be 5 Hz? What am I missing here?